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[57] ABSTRACT 

A computer implemented physical signal analysis method is 
invented. This method includes two essential steps and the 
associated presentation techniques of the results. All the 
steps exist only in a computer: there are no analytic expres- 
sions resulting from the method. The first step is a computer 
implemented Empirical Mode Decomposition to extract a 
collection of Intrinsic Mode Functions (IMF) from 
nonlinear, nonstationary physical signals. The decomposi- 
tion is based on the direct extraction of the energy associated 
with various intrinsic time scales in the physical signal. 
Expressed in the IMF’s, they have well-behaved Hilbert 
Transforms from which instantaneous frequencies can be 
calculated. The second step is the Hilbert Transform. The 
final result is the Hilbert Spectrum. Thus, the invention can 
localize any event on the time as well as the frequency axis. 
The decomposition can also be viewed as an expansion of 
the data in terms of the IMF’s. Then, these IMF’s, based on 
and derived from the data, can serve as the basis of that 
expansion. The local energy and the instantaneous frequency 
derived from the IMF’s through the Hilbert transform give 
a full energy-frequency-time distribution of the data which 
is designated as the Hilbert Spectrum. 
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r 7 30 
otherwise reserves all copyright rights whatsoever. 

BACKGROUND OF THE INVENTION 

1. Technical Field of the Invention 

This invention generally relates to a computer imple- 35 
mented geophysical signal analysis method and apparatus. 
More particularly, this invention relates to a computer imple- 
mented method and apparatus for analyzing nonlinear, non- 
stationary geophysical signals. 

2. Description of Related Art 40 

Analyzing geophysical signals is a difficult problem con- 
fronting many industries. These industries have harnessed 
various computer implemented methods to process data 
taken from geophysical phenomena such as earthquakes, 
ocean waves, tsunamis, ocean surface elevation and wind. 45 
Unfortunately, previous methods have not yielded results 
which are physically meaningful. 

Among the difficulties is that representing geophysical 
processes with geophysical signals may present one or more 5Q 
of the following problems: 

(a) The total data span is too short; 

(b) The data are nonstationary; and 

(c) The data represent nonlinear processes. 

Although problems (a)-(c) are separate issues, the first 55 
two problems are related because a data section shorter than 
the longest time scale of a stationary process can appear to 
be nonstationary. Because many geophysical events are 
transient, the data representative of those events are nonsta- 
tionary. For example, a transient event such as an earthquake 60 
will produce nonstationary data when measured. 
Nevertheless, the nonstationary character of such data is 
ignored or the effects assumed to be negligible. This 
assumption may lead to inaccurate results and incorrect 
interpretation of the underlying physics as explained below. 65 

A variety of techniques have been applied to nonlinear, 
nonstationary geophysical signals. For example, many com- 
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puter implemented methods apply Fourier spectral analysis 
to examine the energy-frequency distribution of such sig- 
nals. 

Although the Fourier transform that is applied by these 
computer implemented methods is valid under extremely 
general conditions, there are some crucial restrictions: the 
system must be linear, and the data must be strictly periodic 
or stationary. If these conditions are not met, then the 
resulting spectrum will not make sense physically. 

A common technique for meeting the linearity condition 
is to approximate the geophysical phenomena with at least 
one linear system. Although linear approximation is an 
adequate solution for some applications, many physical 
phenomena are highly nonlinear and do not admit a reason- 
ably accurate linear approximation. 

Furthermore, imperfect probes/sensors and numerical 
schemes may contaminate data representative of the phe- 
nomenon. For example, the interactions of imperfect probes 
with a perfect linear system can make the final data nonlin- 
ear. 

Many recorded geophysical signals are of finite duration, 
nonstationary, and nonlinear because they are derived from 
geophysical processes that are nonlinear either intrinsically 
or through interactions with imperfect probes or numerical 
schemes. Under these conditions, computer implemented 
methods which apply Fourier spectral analysis are of limited 
use. For lack of alternatives, however, such methods still 
apply Fourier spectral analysis to process such data. 

In summary, the indiscriminate use of Fourier spectral 
analysis in these methods and the adoption of the stationary 
and linear assumptions may give misleading results some of 
which are described below. 

First, the Fourier spectrum defines uniform harmonic 
components globally. Therefore, the Fourier spectrum needs 
many additional harmonic components to simulate nonsta- 
tionary data that are nonuniform globally. As a result, energy 
is spread over a wide frequency range. For example, using 
a delta function to represent the flash of light from a 
lightning bolt will give a phase-locked wide white Fourier 
spectrum. Here, many Fourier components are added to 
simulate the nonstationary nature of the data in the time 
domain, but their existence diverts energy to a much wider 
frequency domain. Constrained by the conservation of 
energy principle, these spurious harmonics and the wide 
frequency spectrum cannot faithfully represent the true 
energy density of the lighting in the frequency and time 
space. 

More seriously, the Fourier representation also requires 
the existence of negative light intensity so that the compo- 
nents can cancel out one another to give the final delta 
function. Thus, the Fourier components might make math- 
ematical sense, but they often do not make physical sense 
when applied. 

Although no physical process can be represented exactly 
by a delta function, some geophysical data such as the near 
field strong earthquake energy signals are of extremely short 
duration. Such earthquake energy signals almost approach a 
delta function, and they always give artificially wide Fourier 
spectra. 

Second, Fourier spectral analysis uses a linear superpo- 
sition of trigonometric functions to represent the data. 
Therefore, additional harmonic components are required to 
simulate deformed wave profiles. Such deformations, as will 
be shown later, are the direct consequence of nonlinear 
effects. Whenever the form of the data deviates from a pure 
sine or cosine function, the Fourier spectrum will contain 
harmonics. Furthermore, both nonstationarity and nonlin- 
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earity can induce spurious harmonic components that cause 
unwanted energy spreading and artificial frequency smear- 
ing in the Fourier spectrum. The consequence is incorrect 
interpretation of physical phenomenon due to the misleading 
energy-frequency distribution for nonlinear and nonstation- 
ary data representing the physical phenomenon. 

Below are several types of nonlinear, nonstationary geo- 
physical signals which are very difficult to analyze with 
traditional computer implemented techniques. The invention 
presented herein, however, is particularly well-suited to 
processing such geophysical signals with a computer. 

Earthquake Signals 

Earthquakes are typically recorded by seismometers such 
as the seismometer 400 which may be implemented with the 
Ranger seismometer manufactured by kinemetrics Model 
WR-1 Wide-Band shown in FIG. 21 which records ground 
accelerations to produce a signal representative of the earth- 
quake. 

Fortunately, all earthquakes are transient lasting only a 
few tenths to a few seconds at most; consequently, earth- 
quake signals are nonstationary. Most earthquake signals are 
still processed with various computer implemented methods 
that apply algorithms based on Fourier analysis (Hu, et al. 
Earthquake Engineering, Chapman & Hall, London, 1996). 
Such earthquake signals are processed to better understand, 
for example, crust structure geophysics, near field earth- 
quakes and site specific ground motions. 

Crust structure geophysics is a term for the geophysical 
structure of the earth which includes the crust and inner core. 
Due to the different geophysical properties of the local crust 
material, the earthquake signal can be used to determine the 
mode of earthquake wave propagation, their dispersion 
characteristics, and the free oscillations. These properties 
can be used to infer the structure of the crust, and the elastic 
properties and density of the crust medium through which 
the wave propagated. 

Most seismologists are interested in the earthquake sig- 
nals to infer the geophysical properties of the earth as 
explained above. Earthquake engineers, however, are inter- 
ested in the destructive power of the earthquakes. Therefore, 
the seismologists prefer sampling the earthquake signal from 
a long distance, up to thousand of miles say, to infer the 
geophysical properties of the crust along the path of wave 
propagation. On the other hand, earthquake engineers are 
most interested in the immediate neighborhood of the earth- 
quake epicenter (near field earthquakes), within a few kilo- 
meters say, where the destruction would be the most severe. 

For any given earthquake, the ground response is site 
specific and depends on the following factors: (1) nature of 
the earthquake (whether it is a shear or a thrust), (2) the 
propagation path, (3) the local ground geo-engineering prop- 
erties (whether it is rock or sediment), and (4) the local 
topographic geometry (e.g., whether in a valley or on the top 
of a hill). These factors influence the severity of the ground 
motion from a given earthquake at specific locations. 

Conventional methods, however, cannot reveal detailed 
information in the dispersion properties, the wave form 
deformation, and the energy-frequency distribution of earth- 
quakes because the data representing the earthquake is 
nonlinear and nonstationary. Revealing this information is 
necessary to correctly understand crust structure geophysics 
and to accurately deduce site specific ground motions. 

Furthermore, most near field strong earthquake ground 
motions are nonstationary because of their extremely short 
duration. Seismic records representing such earthquakes 
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always give artificially wide Fourier spectra because of this 
nonstationarity. This wide frequency distribution will dilute 
the energy content everywhere on the frequency axis and 
distort the true energy-frequency distribution. The result is 
5 that the energy density at critical resonant frequencies for 
specific building structures will be underestimated. The rule 
of thumb for the resonant frequency is given as 1/(0. IN) 
cycle per second, where N is the number of the stories of the 
building. Therefore, for a ten-story building, the resonant 
10 frequency is near 1 Hz. For high-rises, the frequency will be 
even lower. 

For example, the seismic record of the El Centro earth- 
quake is shown in FIG. 10(a). The corresponding Fourier 
spectrum for the El Centro earthquake is illustrated in FIGS. 
15 12(a) and 12(b). The Fourier spectrum shows a conspicuous 
lack of energy in the low frequency range (less than 1 Hz), 
as shown in the detailed spectra covering only 0 to 5 Hz 
frequency range in FIG. 12(b). The low frequency range, 
which is critical to analyzing the effects on the high rise 
20 structures, is severely under-represented. 

This peakiness of energy located at a very narrow fre- 
quency range can cause resonance oscillation of buildings, 
and has been observed at Mexico City during the 1985 

earthquake to have caused great destruction. 

25 

The oscillation could not be duplicated by any linear 
model. Therefore, computer implemented Fourier analysis 
methods are not appropriate for analyzing earthquake sig- 
nals. 

30 Although earthquake wave motion has been approxi- 
mated as a linear process, the strong ground motion may not 
be linear. This nonlinearity is revealed by the distorted wave 
forms in the data (Newmark and Rosenblueth, Fundamen- 
tals of Earthquake Engineering, Prentice Hall, Englewood 
35 Cliffs, N.J, 1971). The rich higher harmonics, a typical 
consequence of Fourier analysis when applied to nonline arly 
deformed signals, have also contributed to clouding the 
characteristics of the real signal and the dynamics of the 
underlying physical systems. 

40 Furthermore, to apply Fourier analysis, some kind of 
window should be used to eliminate the end effects and to 
facilitate the computer implemented computations. For a 
truly nonstationary process, there is no time scale to guide 
the choice of the window size. A long time window (eg. 10 
45 seconds) is necessary for resolving low frequencies (near 0.1 
Hz) that are critical for analyzing the resonant effects on 
modern high-rise structures. Too long a window, however, 
will be plagued by nonstationarity. In other words, within a 
short window the data may look stationary locally. But the 
50 short window will have very poor frequency resolution, for 
the frequency resolution in Hz is given by 1/T, where T is the 
window length. Such difficulties have never been fully 
resolved by applying Fourier spectral analysis to earthquake 
signals. 

55 As explained above, both nonstationarity and nonlinearity 
can induce artificial frequency smearing and obscure the true 
energy density picture and the underlying physics of the 
system being studied. 

6Q Water Wave Signals 

The dynamics of ocean waves are measured from ocean 
sensors located at field stations such as the ocean wave 
sensor 410 which may be implemented by using the NDBC 
3m discus ocean wave sensor illustrated in FIG. 22 which 
65 records an ocean wave signal representing ocean surface 
elevation as a function of time. As example of an ocean wave 
signal is shown in FIG. 13(a). Ocean waves are studied for 
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ship design, ship routing, coastal and off-shore structure 
design, harbor operations, and even weather forecasting. 

Ocean wave signals are typically random and nearly 
nonstationary. In the past, ocean wave signals were analyzed 
by applying computer implemented Fourier analysis. In fact, 5 
the studies of the wave spectra from Fourier analysis have 
been a main subject of ocean wave research (see, for 
example, Huang, et al., 1990a, Wave Spectra, The Sea, 9, 
197-237). 

Traditional computer implemented analysis methods, 10 
however, are not well suited to studying ocean wave signals 
because ocean waves are typically nonlinear and nonstation- 
ary. For example, the corresponding Wavelet spectrum for 
ocean wave signals which is generated by applying a Morlet 
Wavelet transform to the ocean wave signals of FIG. 13 (a) 15 
is shown in FIG. 14(b). The Wavelet spectrum gives a nearly 
continuous distribution, and wide spread of energy consist- 
ing primarily of harmonics in the frequency axis. This 
energy spread is due to the nonlinear and nonstationary 
character of ocean waves. This energy spread also contrib- 20 
utes to the difficulty in analyzing the results of traditional 
computer implemented techniques applying the Wavelet 
transform. 

Water wave signals detected from mechanically generated 25 
water waves by a wave sensor 420 as shown in FIG. 23 have 
been studied to analyze nonlinear water wave evolution 
processes (eg. Huang, et al., The Mechanism for Frequency 
Downshift in Nonlinear Wave Evolution, Advances in 
Applied Mechanics, Vol. 32, pp. 59-117 1996). 30 

Due to weak nonlinear interactions, the frequency of the 
water waves will downshift as they propagate, a process 
necessary for the waves to become longer and grow higher 
under the wind. 

In the narrow-band wave field, the downshift has been 35 
shown as the consequence of the Benjamin-Fier instability 
(Benjamin and Fier, The Disintegration of Wave trains on 
Deep Water, Part I, Theory, J. Fluid Mech., 27, 417-430, 
1967). Although water wave evolution is generally assumed 
to be gradual and continuous, several authors have theorized 40 
that the evolution is not continuous and gradual, but local 
and discrete. 


The resolution power of previous data analysis 
techniques, however, has rendered proof of this theory 
nearly impossible. As explained above, known computer 45 
implemented data analysis techniques which apply Fourier 
analysis are incapable of accurately interpreting nonlinear, 
nonstationary processes such as water wave propagation and 
evolution. 


Tsunami Signals 


Tsunamis are detected with tidal gauges such as the tidal 
gauge 430 shown in FIG. 24 which record water elevation 
as a function of time. An example of a tsunami signal, which 55 
necessarily includes a tidal signal, is shown in FIG. 1 6(a). 

Although tidal signals are generally stationary, tsunami 
waves are transient, nonlinear and nonstationary. Tidal 
gauges necessarily measure both the tide and the tsunami. 
The combined signal, therefore, is nonstationary and non- go 
linear. 

Filtering cannot remove the tsunami signal cleanly 
because the tsunami signals and the tidal signals usually 
have many harmonic components in the same frequency 
range. Therefore, tsunami signals and combined tsunami- 65 
tidal signals lack an effective computer implemented data 
analysis method which is able to handle the nonlinear and 


nonstationary character of the data representative of these 
geophysical processes. 

Ocean Altitude and Ocean Circulation 

Satellite altimetry is a powerful technique for large scale 
ocean circulation studies (Huang, et al. 1978, “Ocean Sur- 
face Measurement Using Elevation From GEOS-3 
Altimeter”, J. Geophys. Res., 83, 4,673^1,682; Robinson, et 
al., 1983, “A Study of the Variability of Ocean Currents in 
the Northwestern Atlantic Using Satellite Altimetry”, J. 
Phys. Oceanogr., 13, 565-585). An orbital satellite system 
440 as shown in FIG. 25 can produce extremely accurate 
data representing the altitude of the ocean surface. 

The accepted view of the equatorial dynamics is the 
propagation of Kelvin waves forced by variable wind stress 
(Byod, 1980, “The Nonlinear Equatorial Kelvin Waves”, J. 
Phys. Oceanogr., 10, 1-11 and Zheng, et al., 1995, “Obser- 
vation of Equatorially Trapped Waves in the Pacific Using 
Geosat Altimeter Data”, Deep-Sea Res., (in press). In this 
model, the wave propagation will leave a surface elevation 
signature of the order of 10 cm, which can be measured by 
the satellite altimeter 440 . 

Because of the importance of the equatorial region in 
determining the global climate pattern, altimeter data have 
been used extensively to study the dynamics of this region 
(Miller, et al., 1988, “GEOSAT Altimeter Observation of 
Kelvin Waves and the 1986-1987 El Nino” Science, 239, 
52-54; Miller, et al., 1990, “Large-Scale Meridional Trans- 
port in the Tropic Pacific Ocean During the 1986-87 El Ni 
no from GEOSAT”, J. Geophys. Res. 95, 17,905-17,919. ; 
Zheng, et al., 1994, “The Effects of Shear Flow on Propa- 
gation of Rossby Waves in the Equatorial Oceans”, J. Phys. 
Oceanogr., 24, 1680-1686 and Zheng, et al., 1995, “Obser- 
vation of Equatorially Trapped Waves in the Pacific Using 
Geosat Altimeter Data”, Deep-Sea Res., (in press)). A typi- 
cal time series on the Equator sea surface elevation data at 
174° East longitude is shown in FIG. 19(a). 

Limited by the data length and complicated by ocean 
dynamics, all the past investigators have faced serious 
problems in processing this nonstationary altimeter data. 
Therefore, weather forecasting which accurately accounts 
for ocean effects has been impossible with traditional com- 
puter implemented data analysis methods. 

According to the above background, the state of the art 
does not provide a useful computer implemented tool for 
analyzing nonlinear, nonstationary data such as geophysical 
signals. 


SUMMARY OF THE INVENTION 

The invention employs a computer implemented Empiri- 
cal Mode Decomposition method which decomposes physi- 
cal signals representative of a physical phenomenon into 
components. These components are designated as Intrinsic 
Mode Functions (IMFs) and are indicative of intrinsic 
oscillatory modes in the physical phenomenon. 

Contrary to almost all the previous methods, this new 
computer implemented method is intuitive, direct, a 
posteriori, and adaptive, with the basis of the decomposition 
based on and derived from the physical signal. The bases so 
derived have no close analytic expressions, and they can 
only be numerically approximated in a specially pro- 
grammed computer by utilizing the inventive methods dis- 
closed herein. 

More specifically, the general method of the invention 
includes two main components or steps to analyze the 
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physical signal without suffering the problems associated 
with computer implemented Fourier analysis, namely inac- 
curate interpretation of the underlying physics caused in part 
by energy spreading and frequency smearing in the Fourier 
spectrum. The first step is to process the data with the 5 
Empirical Mode Decomposition (EMD) method, with which 
the data are decomposed into a number of Intrinsic Mode 
Function (IMF) components. In this way, the signal will be 
expanded by using a basis that is adaptively derived from the 
signal itself. 10 

The second step of the general method of the present 
invention is to apply the Hilbert Transform to the decom- 
posed IMF’s and construct an energy-frequency-time 
distribution, designated as the Hilbert Spectrum, from which 
occurrence of physical events at corresponding times (time 15 
localities) will be preserved. There is also no close analytic 
form for the Hilbert Spectrum. As explained below, the 
invention avoids this problem by storing numerical approxi- 
mations in the specially programmed computer by utilizing 
the inventive method. 20 

The invention also utilizes instantaneous frequency and 
energy to analyze the physical phenomenon rather than the 
global frequency and energy utilized by computer imple- 
mented Fourier spectral analysis. 

Furthermore, a computer implementing the invention, 
e.g., via executing a program in software, to decompose 
physical signals into intrinsic mode functions with EMD and 
generate a Hilbert spectrum is also disclosed. Because of the 
lack of close form analytic expression of either the basis 3Q 
functions and the final Hilbert spectrum; computer imple- 
mentation of the inventive methods is an important part of 
the overall method. 

Still further, the invention may take the form of an article 
of manufacture. More specifically, the article of manufacture 35 
is a computer-usable medium, including a computer- 
readable program code embodied therein wherein the 
computer-readable code causes a computer to execute the 
inventive method. 

Further scope of applicability of the present invention will 40 
become apparent from the detailed description given here- 
inafter. However, it should be understood that the detailed 
description and specific examples, while indicating pre- 
ferred embodiments of the invention, are given by way of 
illustration only, since various changes and modifications 45 
within the spirit and scope of the invention will become 
apparent to those skilled in the art from this detailed descrip- 
tion. Furthermore, all the mathematic expressions are used 
as a short hand to express the inventive ideas clearly and are 
not limitative of the claimed invention. 50 

BRIEF DESCRIPTION OF DRAWINGS 

The present invention will become more fully understood 
from the detailed description given hereinbelow and the 
accompanying drawings which are given by way of illus- 
tration only, and thus are not limitative of the present 
invention, and wherein: 

FIG. 1(a) is a high-level flowchart describing the overall 
inventive method which may be implemented on the com- 6Q 
puter system shown in FIG. 2; 

FIG. 1(b) is a high-level flowchart describing the Sifting 
Process which may be implemented on the computer system 
shown in FIG. 2; 

FIG. 1(c) is a continuation of the high-level flowchart in 65 
FIG. 1(b) describing the Sifting Process which may be 
implemented on the computer system shown in FIG. 2; 
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FIG. 2 is a high-level block diagram of a computer system 
which may be programmed with the inventive method with 
the result being a special purpose computer; 

FIG. 3(a) shows wind speed data in the form of a graph 
plotting wind speed as a function of time for explaining the 
computer implemented Empirical Mode Decomposition 
method of the invention; 

FIG. 3(b) is a graph illustrating the upper envelope, lower 
envelope, envelope mean and original wind speed data 
which are utilized to explain the computer implemented 
Empirical Mode Decomposition method of the invention; 

FIGS. 3(c)-(c) are graphs of the first, second and third 
component signals hi, hll, hl2, respectively which are 
generated by the Sifting Process of the invention; 

FIG. 3(f) is a graph of the first intrinsic mode function 
component which is generated by the Sifting Process of the 
invention; 

FIG. 3(g) is a graph of data with intermittency for 
illustrating an optional intermittency test of the invention; 

FIGS. 3(/z) — (/*) are graphs of the first, second, and third 
intrinsic mode functions when the Sifting Process is applied 
to the data of FIG. 3(g) without applying the intermittency 
test option; 

FIGS. 3 (k)-(m) are graphs of the first, second, and third 
intrinsic mode functions when the Sifting Process is applied 
to the data of FIG. 3(g) which applies the intermittency test 
option. 

FIG. 4 (a) is a graph of a wind speed signal which is for 
explaining the computer implemented Empirical Mode 
Decomposition method of the invention; 

FIGS. 4 (b)-(k) show the wind speed signal and the nine 
intrinsic mode functions which are extracted from the wind 
speed signal by the EMD method of the invention; 

FIGS. 5(a)-(z) are a series of graphs illustrating the 
successive reconstruction of the original wind speed data 
from the intrinsic mode functions; 

FIG. 5(j) illustrates the final difference between the origi- 
nal wind speed data and the reconstructed time series from 
the IMFs. 

FIG. 6(a) is the Hilbert Spectrum generated by the inven- 
tion from the wind speed data of FIG. 4 (a); 

FIG. 6(b) is the conventional Morlet Wavelet spectrum 
generated from the wind speed data of FIG. 4 (a); 

FIG. 6(c) shows the Hilbert Spectrum of FIG. 6(a) after 
smoothing by a 15x15 weighted Gaussian smoothing filter; 

FIG. 7 is a comparison of the marginal Hilbert spectrum 
(solid line) and the fourier spectrum (dotted line) which 
were generated from the wind speed signal of FIG. 4 (a); 

FIG. 8(a) is a graph illustrating the Degree of Stationarity 
and Degree of Statistical Stationarity which were generated 
from the wind speed signal of FIG. 4(a) with time averages 
of 10, 50, 100 and 300; 

FIGS. 8(b) and (c) are sections of the wind speed data that 
was used by the invention to produce the Degree of Sta- 
tionarity shown in FIG. 8(a); 

FIG. 9(a) shows the profile of a Stokes wave in deep water 
which may be processed by the invention; 

FIG. 9(b) shows the IMF generated by the invention from 
the Stokes wave shown in FIG. 9(a); 

FIG. 9(c) is a conventional Morlet Wavelet spectrum of 
the Stokes wave shown in FIG. 9(a); 

FIG. 9(d) is the Hilbert Spectrum generated by the inven- 
tion from the Stokes wave shown in FIG. 9(a); 
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FIG. 10(a) is an earthquake signal representative of the El 
Centro earthquake; 

FIGS. 10(6)-(6) show the ten intrinsic mode functions 
which are extracted from the El Centro earthquake signal by 
the EMD method of the invention; 

FIG. 11(a) is the Hilbert Spectrum generated by the 
invention from the El Centro earthquake signal shown in 
FIG. 10(a); 

FIG. 11(6) is a conventional Morlet Wavelet spectrum of 
the El Centro earthquake signal shown in FIG. 10(a); 

FIG. 11(c) shows the Hilbert Spectrum of FIG. 11(a) after 
smoothing by a 15x15 weighted Gaussian smoothing filter; 

FIGS. \2(a)-(d) compare the marginal Hilbert spectrum 
with the fourier spectrum of the El Centro earthquake signal 
with FIG. 12(a) showing the fourier spectrum with full 
frequency range, FIG. 12 (6) showing the fourier spectrum 
for the detailed 5 Hz range, FIG. 12(c) showing the marginal 
spectrum with full frequency range, FIG. 12(d) showing the 
marginal spectrum for the detailed 5 Hz range; 

FIG. 13(a) is an ocean wave signal taken from a tidal 
station; 

FIGS. 13(6)-0 show the eight intrinsic mode functions 
which are extracted from the ocean wave signal of FIG. 
13(a) by the EMD method of the invention; 

FIG. 14(a) is the Hilbert Spectrum generated by the 
invention from the ocean wave signal shown in FIG. 13(a) 
after smoothing by a 9x9 weighted Gaussian smoothing 
filter; 

FIG. 14(b) is a conventional Morlet Wavelet spectrum of 
the ocean wave signal shown in FIG. 13 (a); 

FIG. 14(c) shows the Hilbert Spectrum from a 5-minute 
section of the ocean wave signal plotted together with the 
ocean wave signal in an arbitrary scale; 

FIG. 14(d) shows a conventional Wavelet spectrum the 
same 5-minute section of the ocean wave signal from FIG. 
14(c) plotted together with the ocean wave signal in an 
arbitrary scale; 

FIG. 14(c) compares the Wavelet and Hilbert spectra of 
the ocean wave signal; 

FIG. 140 compares the fourier spectrum (bottom line), 
marginal Hilbert spectrum (middle line) and Wavelet spec- 
trum (top line) generated from the ocean wave signal; 

FIG. 15(a) is a water wave signal representing a water 
surface elevation changes from mechanically generated 
water waves; 

FIG. 15(6) shows the unwrapped phase functions gener- 
ated from the water wave signal of FIG. 15(a); 

FIG. 15(c) shows the phase change relative to the values 
at station 1 which senses the water wave signal shown in 
FIG. 15(a); 

FIGS. 15(d)-(g) are phase amplitude diagrams for sta- 
tions 2,4,5 and 6, respectively each of which senses the 
water wave signal shown in FIG. 15(a) at a different 
location; 

FIG. 15(6) is a two way comparison of a water wave 
profile change at the location of a phase jump; 

FIG. 15(/) shows a Hilbert Spectrum of a water wave 
signal for the purpose of illustrating water wave evolution 
analysis according to the invention; 

FIG. 150 shows a marginal spectrum of a water wave 
signal according to the invention; 

FIG. 15(6) is a conventional fourier spectrum of the same 
water wave signal processed in FIG. 150; 
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FIG. 150 a conventional Morlet Wavelet spectrum of the 
same water wave signal processed in FIG. 150; 

FIG. 15(m) is a water wave signal plotted together with its 
Hilbert spectrum which is generated by the invention; 

5 FIG. 150) is a water wave signal plotted together with a 
conventional Wavelet Spectrum; 

FIG. 16(a) is a combined tidal signal and tsunami signal 
collected inside the Kahului Harbor in Maui Hawaii from 
Oct. 4 to Oct. 9, 1994; 

10 

FIGS. 16(6)-(z) show the eight intrinsic mode functions 
which are extracted from the combined tidal signal and 
tsunami signal of FIG. 16(a) by the EMD method of the 
invention; 

15 FIG. 17 shows the separation of the tidal signal and the 
tsunami signal by the IMF components generated by the 
invention; 

FIG. 18(a) shows the corresponding Hilbert Spectrum for 
the combined tidal signal and tsunami signal of FIG. 16(a); 
20 FIG. 18(6) shows the corresponding Hilbert Spectrum for 
the tsunami signal component; 

FIG. 18(c) shows the corresponding marginal Hilbert 
spectrum for the tidal signal component; 

FIG. 18(0 shows the marginal spectrum of the tidal signal 
25 component; 

FIG. 19(a) shows ocean surface elevation signal taken by 
an orbiting satellite altimeter at 174E on the equator on Apr. 
1, 1985; 

FIG. 19(6) shows the conventional Morlet Wavelet Spec- 
30 trum generated from the ocean surface elevation signal 
shown in FIG. 19(a); 

FIGS. 19(c)-(g) show the five intrinsic mode functions 
which are extracted from the ocean surface elevation signal 
shown in FIG. 19(a) by the EMD method of the invention; 

FIG. 19(6) shows the Hilbert Spectrum generated from 
the ocean surface elevation signal shown in FIG. 19(a) by 
the invention; 

FIG. 20(a) is a special Hilbert Spectrum with only the 
40 energy containing IMFs (c2 to c6) generated from a wind 
speed signal; 

FIG. 20(6) shows the instantaneous energy density of the 
fluctuating components of the wind speed signal (dotted 
line) and the overall energy density with the trend (solid 
45 line) which are generated by the invention; 

FIG. 21 shows a typical seismometer detecting an earth- 
quake and generating an earthquake signal representative 
thereof; 

FIG. 22 shows an ocean wave sensor in a field station for 
50 detecting an ocean wave and generating an ocean wave 
signal indicative thereof; 

FIG. 23 shows a wave sensor detecting mechanically 
generated waves and generating a wave signal indicative 
thereof; 

55 FIG. 24 shows a tidal gauge detecting tidal fluctuations 
and tsunamis and generating a combined tidal signal and 
tsunami signal indicative thereof; and 

FIG. 25 shows an orbital satellite equipped with an 
altimeter detecting ocean surface elevation and generating 
60 an ocean surface elevation signal indicative thereof. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

Before describing the computer implemented Empirical 
65 Mode Decomposition method in detail, the definition and 
physical meaning of intrinsic mode functions will be dis- 
cussed. 
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Intrinsic Mode Function 

A Intrinsic Mode Function (IMF) is a function that 
satisfies the following two conditions: 

(a) in the whole data set, the number of extrema and the 
number of zero -crossings must either be equal or differ 
at most by one, and 

(b) at any point, the mean value of the envelope defined 
by the local maxima and the envelope defined by the 
local minima is zero. 

The first condition shares some similarity to the tradi- 
tional narrow band requirements for a stationary Gaussian 
process. The second condition is a totally new idea. 
Conceptually, the second condition modifies the classical 
global requirement to a local one. Furthermore, the second 
condition has the desirable result that the instantaneous 
frequency will not have unwanted fluctuations induced by 
asymmetric wave forms. Mathematically, the second condi- 
tion should ideally be The local mean of the data being zero/ 
For nonstationary data, the Tocal mean’ requires a Tocal 
time scale’ to compute the mean, which is impossible to 
define. Fortunately, the local time scale need not be defined 
to fulfil the second condition, as will be discussed below. 

To apply these concepts to physical data, the invention 
utilizes the local mean of the signal envelopes to force the 
local symmetry. The signal envelopes are defined by the 
local maxima and the local minima. This is an approxima- 
tion which avoids the definition of a local averaging time 
scale. With the physical approach and the approximation 
adopted here, the inventive method does not always guar- 
antee a perfect instantaneous frequency under all conditions. 
Nevertheless, it can be shown that, even under the worst 
conditions, the instantaneous frequency so defined is still 
consistent with the physics of the system being studied and 
represents the system being studied much more accurately 
than previous techniques based on Fourier analysis. 

The term “Intrinsic Mode Function” is adopted because it 
represents the oscillation mode embedded in the data. With 
this definition, the IMF in each cycle, defined by the 
zero-crossings, involves only one mode of oscillation. In 
other words, each IMF represents only one group of oscil- 
lation modes or time scales and no riding waves are allowed. 

Before presenting the inventive EMD method for decom- 
posing the data into IMFs, a qualitative assessment of the 
intrinsic oscillatory modes may be roughly determined by 
simply examining the data by eye. From this examination, 
one can immediately identify the different scales directly in 
two ways: the time lapse between the successive alternations 
of local maxima and minima and the time lapse between the 
successive zero-crossings reveals the different scales. The 
interlaced local extrema and zero -crossings give us compli- 
cated data: one undulation is riding on top of another, and 
they, in turn, are riding on still other undulations, and so on. 
Each of these undulations defines a characteristic scale or 
oscillation mode that is intrinsic to the data: hence, the term 
“Intrinsic Mode Function” is adopted. 

To reduce the data into the needed IMFs, the invention 
utilizes a computer implemented Empirical Mode Decom- 
position Method which is described below. 

Empirical Mode Decomposition (EMD): The 
Sifting Process 

First, the Empirical Mode Decomposition method which 
deals with both nonstationary and nonlinear data will be 
discussed. Then, the physical meaning of this decomposition 
will be presented. 

The essence of the EMD method is to identify empirically 
the intrinsic oscillatory modes by their characteristic time 
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scales in the data, and then decompose the data accordingly. 
The decomposition is based on the following assumptions: 

a. the signal has at least two extrema: one maximum and 
one minimum, and 

5 b. the characteristic time scale is defined by the time lapse 
between the extrema. 

In other words, the invention adopts the time lapse 
between successive extrema as the definition of the time 
scale for the intrinsic oscillatory mode because it gives a 
10 much finer resolution of the oscillatory modes and because 
it can be applied to data with non-zero mean (either all 
positive or all negative values, without zero -crossings). A 
systematic way to extract the intrinsic mode functions is the 
computer implemented Empirical Mode Decomposition 
15 method or Sifting Process which is described as follows. 

FIG. 1(a) illustrates the overall inventive method includ- 
ing the Sifting Process in step 120. First, the physical 
activity, process or phenomenon is sensed by an appropriate 
sensor in step 100. Appropriate sensors for detecting the 
20 physical activity and generating a physical signal represen- 
tative thereof are discussed in the practical examples below. 

After sensing in step 100, the analog signal is converted 
to the digital domain suitable for computer processing in the 
A/D conversion step 105. 

25 Next, an optional smoothing step 110 may be applied to 
the physical signal. The optional smoothing step 110 may be 
applied to smooth the signal with, for example, a weighted 
running average to remove excessive noise. 

Thereafter, the Sifting Process is applied in step 120 to sift 
30 the signal with the Empirical Mode Decomposition method 
and thereby extract the intrinsic mode function(s). The 
intrinsic mode functions can then be displayed as shown in 
step 130 and checked for orthogonality in step 135. 

Before continuing with the main flow in FIG. 1(a ) , the 
35 details of the Sifting Process will be explained with refer- 
ence to the high level flowchart in FIGS. 2(a ) , 2(b) and the 
series of graphs showing illustrative results of the Sifting 
Process in FIGS. 3 («)-(/)• 

As shown in FIG. 1(b), the digitized physical signal from 
40 step 105 is first windowed by framing the end points in step 
107. Then, the Sifting Process begins at step 200 by iden- 
tifying local maximum values of the digitized, framed 
physical signal from step 107. FIG. 3(a) shows a typical 
physical signal 10 which, in this example, represents wind 
45 speed spanning a time interval of one second. 

Before construction of the envelope in steps 210 and 230, 
optional intermittency tests ( 201 , 221 ) may be introduced to 
alleviate the alias associated with intermittence in the data 
that can cause mode mixing. 

50 Optional intermittency test 201 checks the distance 
between successive maxima to see if this distance between 
is within a pre- assigned value n times the shortest distance 
between waves. If no, then an intermittency exists and the 
method proceeds to step 202. If yes, then there is no 
55 intermittency and the upper envelope is constructed in step 
210 as further described below. 

Similarly optional intermittency test 221 checks the dis- 
tance between successive minima to see if this distance is 
within a pre -assigned value n times the shortest distance 
60 between waves. If no, then an intermittency exists and the 
method proceeds to step 222. If yes, then there is no 
intermittency and the upper envelope is constructed in step 
230 as further described below. 

An example of such intermittency is given in FIG. 3(g), 
65 in which small scale waves appear only intermittently. By 
strict application of the Sifting Process, the resulting IMFs 
are given in FIGS. 3 (b)-(j), in which two drastically differ- 
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ent time scales are present in the first IMF component as 
shown in FIG. 3(h). This mixing of modes breaks up the 
main wave train by the small intermittent oscillations. 

If intermittency tests (201,222) are employed which uti- 
lize a preassigned value of n-times the shortest distance 
between waves, the resulting IMFs are shown in FIGS. 
3 (k)-(m), in which the modes are clearly and physically 
separated. The effective step to eliminate the mode mixing 
is achieved by treating the local extrema which failed the 
intermittency test as local maxima and minima (steps 202 
and 212), respectively. Therefore, the upper and lower 
envelope will be identical as the original data reference line. 

These intermittency tests (201,221) and the further steps 
(202,222) are optional. By selecting an artificially large n 
value utilized in steps 201 and 221 to test for intermittency, 
the test will be effectively passed. Otherwise, the test can be 
bypassed at the initial selection of the program. 

Then, the method constructs an upper envelope 20 of the 
physical signal 10 in step 210. The upper envelope 20 is 
shown in FIG. 3(b) using a dot-dashed line. This upper 
envelope 20 is preferably constructed with a cubic spline 
that is fitted to the local maxima. 

Next, the local minimum values of the physical signal 10 
are identified in step 220. To complete the envelope, a lower 
envelope 30 is constructed from the local minimum values 
in step 230. The lower envelope 30 is also shown in FIG. 
3(b) using a dot-dash line. 

Like the upper envelope 20, the lower envelope 30 is 
preferably constructed with a cubic spline that is fitted to the 
local minima. 

The upper and lower envelopes 20, 30 should encompass 
all the data within the physical signal 10. From the upper and 
lower envelopes 20, 30, an envelope mean 40 is determined 
in step 240. The envelope mean 40 is the mean value of the 
upper and lower envelopes 20, 30. As can be seen in FIG. 
3(b), the envelope mean 40 bisects the physical signal 10 
quite well. 

Then, the method generates the first component signal h ± 
in step 250 by subtracting the envelope mean 40 from the 
physical signal 10. This computer implemented step may 
also be expressed as: 

( 1 ) 

Where the envelope mean 40 is m 1 and the physical signal 
10isX(t). 

FIG. 3(c) shows the first component signal h ± . Ideally, the 
first component signal h 2 should be an IMF, for the con- 
struction of h 1 described above seems to have made h 1 
satisfy all the requirements of an IMF. In reality, however, a 
gentle hump that resides on a slope region of the data can 
become an extremum when the reference coordinate is 
changed from the original rectangular coordinate to a cur- 
vilinear coordinate. For example, the imperfection of the 
envelopes 20, 30 can be seen by observing the overshoots 
and undershoots at the 4.6 and 4.7 second points in FIG. 
3(b). 

An example of this amplification can be found in the 
gentle hump between the 4.5 and 4.6 second range in the 
data in FIG. 3(a). After the first round of sifting, the gentle 
hump becomes a local maximum at the same time location 
in the first component signal h 2 shown in FIG. 3(c). New 
extrema generated in this way actually recover the proper 
modes lost in the initial examination. Thus, the Sifting 
Process extracts important information from the signal 
which may be overlooked by traditional techniques. In fact, 
the Sifting Process can recover low amplitude riding waves, 
which may appear as gentle humps in the original signal, 
with repeated siftings. 
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Still another complication is that the envelope mean 40 
may be different from the true local mean for nonlinear data. 
Consequently, some asymmetric wave forms can still exist 
no matter how many times the data are sifted. This must be 
5 accepted because, after all, the inventive method is an 
approximation as discussed before. 

Other than these theoretical difficulties, on the practical 
side, serious problems of the spline fitting can occur near the 
ends, where the cubic spline fitting can have large swings. 
10 Left by themselves, the end swings can eventually propagate 
inward and corrupt the whole data span especially in the low 
frequency components. A numerical method has been 
devised to eliminate the end effects details of which will be 
given later. Even with these problems, the Sifting Process 
15 can still extract the essential scales from the data. 

The Sifting Process serves two purposes: to eliminate 
riding waves and to make the wave profiles more symmetric. 
Toward these ends, the Sifting Process has to be repeated. 
Because only the first component signal h 2 has been gener- 
20 ated so far, the decision step 260, which tests successive 
component signals to see if they satisfy the definition of an 
IMF, is bypassed during the first iteration. 

Thus, step 265 is performed which treats the component 
signal as the physical signal in the next iteration. The next 
25 iteration is then performed by executing steps 200-250. In 
step 250, the second component signal h 1:L is generated by 
subtracting the envelope mean from the physical signal (in 
this iteration, the first component signal h 1 is treated as the 
physical signal). In more formal terms: 

30 

^i -w ii = ^ii- (2) 

FIG. 3(d) shows the second component signal h 1± . 
Although the second sifting shows great improvement in the 
signal with respect to the first sifting, there is still a local 
35 maximum below the zero line. After a third sifting, the result 
(third component signal h 12 ) is shown in FIG. 3(d). Now all 
the local maxima are positive, and all the local minima are 
negative, but to ensure this configuration is stable, the 
Sifting Process should be further repeated. In general, the 
40 Sifting Process is repeated at least 3 more times and, in 
general, K times to produce h lk . If no more new extrema are 
generated, then h lk is an IMF. In formal terms: 

K(k-iT m ik= h ik’ ( 3 ) 

45 

The resulting first IMF component is shown in FIG. 3(f) 
after 9 siftings. The first IMF component of the physical 
signal may be designated as such in step 270 and stored in 
step 275 in memory 415: 

50 Cl =h lk , ( 4 ) 

As mentioned above, all these manipulations are carried 
out numerically in a computer 410. There is no explicit close 
form analytic expression for any of the computer imple- 
55 mented steps. 

As described above, the process is indeed like sifting of 
the data by the computer 410 because it separates the finest 
(shortest time scale) local mode from the data first based 
only on the characteristic time scale. The Sifting Process, 
60 however, has two effects: 

a. to eliminate riding waves, and 

b. to ensure the envelopes generated by maxima and 
minima are symmetric. 

While the first condition is necessary for the instantaneous 
65 frequency to be meaningful (as discussed below), the second 
condition is also necessary in case the neighboring wave 
amplitudes have too large a disparity. 
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Unfortunately, the effects of the second condition, when 
carried to the extreme, could obliterate the physically mean- 
ingful amplitude fluctuations. Therefore, the Sifting Process 
should be applied with care, for carrying the process to an 
extreme could make the resulting IMF a pure frequency 
modulated signal of constant amplitude. 

To guarantee that the IMF component retains enough 
physical sense of both amplitude and frequency 
modulations, a stopping criterion is employed to stop the 
generation of the next IMF component. 

This stopping criterion is part of the computer imple- 
mented method and is shown as step 260 in FIG. 1(c). Step 
260 is a decision step that decides whether the stopping 
criteria has been satisfied. The preferred stopping criteria 
determines whether three successive component signals 
satisfy the definition of IMF. If three successive component 
signals all satisfy the definition of IMF, then the Sifting 
Process has arrived at an IMF and should be stopped by 
proceeding to step 270. If not, step 260 starts another 
iteration by proceeding to step 265 as described above. 

Alternatively, another stopping criteria could be used that 
determines whether successive component signals are sub- 
stantially equal. If successive component signals are sub- 
stantially equal, then the Sifting Process has arrived at an 
IMF and should be stopped by proceeding to step 270. If not, 
step 260 starts another iteration by proceeding to step 265 as 
described above. 

Determining whether successive component signals are 
substantially equal in the alternative stopping criteria limits 
the size of the standard deviation, sd, computed from the two 
consecutive sifting results as follows: 

T ( 5 ) 

a<= y| l(*m-i)(0-*u(0)l 2 

2-M (0 
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residual signal as the physical signal in the next iteration. 
Thereafter, the next iteration is performed beginning with 
the execution of step 200 as described above. 

The Sifting Process is analogous to a mechanical sieve, 
except it is implemented here in specially programmed 
computer and applied to any digital data numerically rather 
than mechanically. 

The Sifting Process is repeated for all the subsequent r 7 -’s. 
This iterative procedure may be expressed as: 

r i~ c 2 =r 2 , 


r„-i ~c n =r n . ( 7 ) 

Step 300 stops the Sifting Process by proceeding to stop 
step 310 if the residual signal r n has no more than 2 extrema. 
Otherwise, the method proceeds to step 320. 

In other words, Step 310 stops the Sifting Process if the 
residual signal r„ is monotonically increasing or decreasing. 
This stopping criterion is based on the fact that an IMF 
cannot be extracted from a monotonic function. If r n is not 
monotonically increasing/decreasing, then a next iteration is 
performed beginning with step 320. 

Even for data with zero mean, the final residue still can be 
different from zero. For data with a trend, the final residue 
should be that trend. 

In summary, the Sifting Process decomposes the physical 
signal X(t) into a series of intrinsic mode functions and a 
residue which may be expressed as: 

A (8) 

X(t) = ^ c i + r n . 


Avery rigorous and preferred value for sd is set is desired, 
then a trade-off such as a less rigorous value for sd may be 
used. 

Overall, the first IMF component c 1 should contain the 
finest scale or the shortest period component of the physical 
signal 10. 

Before extracting the next IMF component, a test should 
be conducted to determine if the Sifting Process should stop. 
The stopping criteria is shown in Step 300. Step 300 
determines whether the component signal has more than 2 
extrema. If not, all of the IMF’s have been extracted and the 
Sifting Process is stopped by proceeding to step 310. If so, 
then additional IMF’s may be extracted by continuing the 
process in step 320. 

Step 270 recognizes that an IMF component has been 
successfully generated by the Sifting Process by setting the 
component signal equal to an intrinsic mode function. More 
specifically, step 270 causes the computer 410 to store the 
component signal hlk as an intrinsic mode function in 
memory 415. 

Then, the first IMF is separated from the physical signal 
in step 290 to generate a residual signal. In particular, a 
residual signal is generated by subtracting the intrinsic mode 
function from the physical signal. In formal terms: 

m-ci=r ± . ( 6 ) 

Because the residue, r 1? still includes information of 
longer period components, it is treated as the new physical 
data and subjected to the same Sifting Process as described 
above. Step 320 performs this function by treating the 


In other words, the invention extracts a series of IMFs by 
sifting the physical signal with a computer implemented 
Empirical Mode Decomposition method. This IMF series 
cannot be generated or derived by any analytic method. It 
can only be extracted by the invention through a specially 
programmed computer through the Sifting Process. 

Computer for Implementing Inventive Method 

A computer suitable for programming with the inventive 
method is diagrammatically shown in the block diagram of 
FIG. 2. The computer 410 is preferably part of a computer 
system 400. To allow human interaction with the computer 
410, the computer system includes a keyboard 430 and 
mouse 435. The computer programmed with the inventive 
method is analogous to a mechanical sieve: it separates 
digital data into series of IMF’s according to their time 
scales in a manner analogous to a mechanical sieve which 
separates aggregated sand according to their physical size. 

Because the invention is applied to analyze physical 
signals, the computer system 400 also includes an input 
device 440, sensor 442 and/or probe 444 which are used to 
sample a physical phenomenon and generate physical signal 
representative thereof. More particular examples of such 
inputs (440, 442 and 444) are described in relation to FIGS. 
21-25. 

To output the results of the computer implemented 
method, the computer system 400 also includes a display 
450 such as a cathode ray tube or flat panel display, printer 
460 and output device 470. Each of these outputs (450, 460, 
470) should have the capability to generate color outputs 
because, for example, the Hilbert Spectrum may be in color. 



5,983,162 


17 

Furthermore, the computer system 400 also includes a 
mass storage device 420. The mass storage device 420 may 
be a hard disk, floppy disc, optical disc, etc. The mass 
storage device 420 may be used to store a computer program 
which performs the invention when loaded into the com- 
puter 410. As an alternative, the input device 440 may be a 
network connection or off-line storage which supplies the 
computer program to the computer 410. 

More particularly, the computer program embodiment of 
the invention may be loaded from the mass storage device 
420 into the internal memory 415 of the computer 410. The 
result is that the general purpose computer 410 is trans- 
formed into a special purpose machine that implements the 
invention. 

Even more particularly, each step of inventive method 
will transform at least a portion of the general purpose 
computer 410 into a special purpose computer module 
implementing that step. For example, when the sifting step 
120 is implemented on the computer 410, the result is a 
computer implemented sifting apparatus (sifter) that per- 
forms the sifting functions of sifting step 120 . 

Other embodiments of the invention include firmware 
embodiments and hardware embodiments wherein the 
inventive method is programmed into firmware (such as 
EPROM, PROM or PLA) or wholly constructed with hard- 
ware components. Constructing such firmware and hardware 
embodiments of the invention would be a routine matter to 
one of ordinary skill using known techniques. 

Article of Manufacture 

Still further, the invention disclosed herein may take the 
form of an article of manufacture. More specifically, the 
article of manufacture is a computer-usable medium, includ- 
ing a computer-readable program code embodied therein 
wherein the computer-readable code causes computer 410 to 
execute the inventive method. 

A computer diskette such as disc 480 in FIG. 2 is an 
example of such a computer-usable medium. When the disc 
480 is loaded into the mass storage device 480, the 
computer-readable program code stored therein is trans- 
ferred into the computer 410. In this way, the computer 410 
may be instructed to perform the inventive methods dis- 
closed herein. 

Illustration of Sifting Process 

To further illustrate the Sifting Process, consider FIG. 
4 (a) which shows a physical signal representing wind data 
collected in a laboratory wind-wave tunnel with a high 
frequency response Pitot tube located 10 cm above the mean 
water level. The wind speed was recorded under the condi- 
tion of an initial onset of water waves from a calm surface. 
Clearly, the physical signal is quite complicated with many 
local extrema but no zero -crossings such that the time series 
represents all positive numbers. 

Although the mean can be treated as a zero reference, 
defining it is difficult, for the whole process is transient. This 
example illustrates the advantage of adopting the successive 
extrema for defining the time scale and the difficulties of 
dealing with noristationary data. In fact, a physically mean- 
ingful mean for such data is impossible to define using 
standard methods, the EMD eliminates this difficulty. 

FIG. 4(b) shows the wind speed signal of FIG. 4(a) on a 
different scale for comparison purposes. FIGS. 4 (c)-(k) 
show all the IMFs obtained from repeatedly sifting the wind 
speed signal in FIG. 4(b). The efficiency of the invention is 
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also apparent: the Sifting Process produces a total of 9 
intrinsic mode function components while the Fourier trans- 
form needs components which number as many as half of 
the total number of points in a given window to represent the 
5 wind data with the same accuracy. 

The separation of the wind speed data into locally non- 
overlapping time scale components is clear from FIGS. 
4 (c)-(k). In some components, such as c 2 and c 3 , the signals 
are intermittent, then the neighboring components might 
10 include oscillations of the same scale, but signals of the 
same time scale would never occur at the same locations in 
two different IMF components. 

The components of the EMD are usually physical, for the 
characteristic scales are physically meaningful. 

To confirm the validity and completeness of the Sifting 
Process, the wind speed signal can be reconstructed from the 
IMF components. FIGS. 5(a)-(z) show this reconstruction 
process starting from the longest period IMF to the shortest 
20 period IMF in sequence. For example, FIG. 5(a) shows the 
wind speed signal and the longest period component, C 9 , 
which is actually the residue trend, not an IMF. 

By itself, the fitting of the trend is quite impressive, and 
it is very physical: the gradual decrease of the mean wind 
25 speed indicates the lack of drag from the calm surface 
initially and the increasing of drag after the generation of 
wind waves. As the mean wind speed deceases, the ampli- 
tude of the fluctuation increases, another indication of wind- 
wave interactions. 

30 By adding the next longest period component, c 8 , the 
trend of the sum, c 9 +c 8 , takes a remarkable turn, and the 
fitting to the wind speed signal looks greatly improved as 
shown in FIG. 5(b). Successively adding more components 
with increasing frequency results in the series of FIGS. 
35 5(c)-(z). The gradual change from the mono tonic trend to the 
final reconstruction is illustrative by itself. By the time the 
sum of IMF components reaches C 3 in FIG. 5(g), essentially 
all the energy contained in the wind speed signal is recov- 
ered. The components with the highest frequencies add little 
40 more energy, but they make the data look more complicated. 
In fact, the highest frequency component is probably not 
physical, for the digitizing rate of the Pitot tube is too slow 
to capture the high frequency variations. As a result, the data 
are jagged artificially by the digitizing steps at this fre- 
45 quency. The difference between the original data and the 
re -constituted set from the IMF’s is given in FIG. 5(f). The 
magnitude of the difference is 10 -15 , which is the limit of the 
computer 410. 

50 The Hilbert Spectrum 

Having obtained the Intrinsic Mode Function 
components, the next step in the computer implemented 
method is to apply the Hilbert Transform to each component 
and generate the Hilbert Spectrum as shown in step 140 in 
55 fig. 1(a). 

A brief tutorial on the Hilbert transform with emphasis on 
its physical interpretation can be found in Bendat and 
Piersol, 1986, “Random Data: Analysis and Measurement 
Procedures ”, 2nd Ed., John Wiley & Sons, New York, N.Y. 
60 Essentially, the Hilbert 



65 

transform is the convolution of X(t) with 1/t. This convo- 
lution emphasizes the local properties of X(t). In more 
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formal terms, given a time series X(t), the Hilbert Transform 
Y(t) can be expressed as where P indicates the Cauchy 
principal value. 

With this definition, X(t) and Y(t) form a complex con- 
jugate pair. This complex conjugate pair Z(t) may be 5 
expressed as: 

Z(t)=X(t)+iY(t)=a(t)e i8(l \ (10) 

in which 
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After performing the Hilbert transform to each IMF 
component, we can express the time series data 



In Equation (13), the residue, r n , is purposefully omitted, 
for it is either a monotonic function, or a constant. Although 
the Hilbert transform can treat the monotonic trend as part 25 
of a longer oscillation, the energy involved in the residual 
trend could be overpowering. In consideration of the uncer- 
tainty of the longer trend, and in the interest of information 
contained in the other low energy and higher frequency 
components, the final non-IMF component should be left 30 
out. It, however, could be included, if physical consider- 
ations justify its inclusion. 

Note that Equation (13) gives both amplitude and fre- 
quency of each component as functions of time. It should be 
pointed out that no analytical method can generate the 35 
expression in Equation (13). Instead, all the components 
may be extracted only by a specially programmed computer 
applying the inventive Sifting Process and the Hilbert trans- 
form. The variable amplitude and frequency have not only 
greatly improved the efficiency of the expansion, but also 
enabled the expansion to accommodate nonstationary data. 40 
With IMF expansion, the amplitude and the frequency 
modulations are also clearly separated. 

Equation (13) also enables the computer implemented 
method to represent the amplitude and frequency as func- 
tions of time in a three-dimensional plot, in which the 45 
amplitude can be contoured on the frequency-time plane. 
This frequency-time distribution of the amplitude is desig- 
nated as the Hilbert Amplitude Spectrum, H(co, t), or simply 
Hilbert Spectrum. Thus we have: 

50 

n (14) 

H(io, t) = 

j= 1 


In which H(co, t) is the Hilbert spectrum of the frequency ( 00 ) 55 
and time (t) and a 7 -(t) is the j-th component of the IMF. In the 
presentation, the amplitude (with or without smoothing) can 
be expressed in color maps, black-grey maps, or contour 
maps. Color maps, however, greatly increase the operator’s 
ability to fully analyze the spectrum. In some cases, a color 60 
map will permit the operator to discern relationships and 
trends that would not be apparent in black-grey maps 
thereby making a color display a necessary component in 
some cases. 

If amplitude squared is more desirable to represent energy 65 
density, then the squared values of amplitude can be sub- 
stituted to produce a Hilbert Energy Spectrum just as well. 
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Various forms of Hilbert spectra presentations can be 
generated by the computer in the display step 190: both color 
coded maps and contour maps may be employed to present 
the Hilbert spectra with or without smoothing. The Hilbert 
Spectrum in the color map format for the wind data is shown 
in FIG. 6(a). The Hilbert spectrum in FIG. 6(a) gives a very 
different appearance when compared with the corresponding 
Wavelet spectrum shown in FIG. 6(b). While the Hilbert 
Spectrum in FIG. 6(a) appears only in the skeleton form 
with emphasis on the frequency variations of each IMF, the 
Wavelet analysis result gives a smoothed energy contour 
map with a rich distribution of higher harmonics. 

If a more continuous form of the Hilbert Spectrum is 
preferred, a smoothing method can be optionally applied in 
step 155. The first type of a smoothing method which may 
be used in the invention is a weighted spatial filter which 
averages over a range of cells. For example, a 15 by 15 
weighted Gaussian filter may be employed in step 155 as is 
known in the art to smooth this data. FIG. 6(c) shows the 
result of applying the 15 by 15 weighted Gaussian filter. 

Although smoothing step 155 degrades both frequency 
and time resolutions, the energy density and its trends of 
evolution as functions of frequency and time are easier to 
identify. In general, if more quantitative results are desired, 
the original skeleton presentation is better. If more qualita- 
tive results are desired, the smoothed presentation is better. 
As a guide, the first look of the data is better in the smoothed 
format. 

The alternative of the spatial smoothing is to select a 
lower frequency resolution and leave the time axis undis- 
turbed. The advantages of this approach are the preservation 
of events’ locations and a more continuous frequency varia- 
tion. Furthermore, a lower frequency resolution saves com- 
putation time for the computer implemented method. 

To optimize such computation time, the optimal fre- 
quency resolution in the Hilbert spectrum can be computed 
as follows. Let the total data length be T, and the digitizing 
rate of the sensor be At. Then, the lowest frequency that can 
be extracted from the data is 1/T Hz, which is also the limit 
of frequency resolution for the data. The highest frequency 
that can be extracted from the data is l/(n At) Hz, in which 
n represents the minimum number of At needed to define the 
frequency accurately. 

Because the Hilbert transform defines instantaneous fre- 
quency by differentiation, more data points are needed to 
define an oscillation. The absolute minimum number of data 
points is five for a whole sine wave. Although a whole sine 
wave is not needed to define its frequency, many points 
within any part of the wave are needed to find a stable 
derivative. Therefore, the maximum number of the fre- 
quency cells, N, of the Hilbert spectrum should be 

±_ (15) 

N= nto = ZL 

1 

T 

In order to make the derivative stable, the data is averaged 
over three adjacent cell values for the final presentation. 

To illustrate, consider the wind data of FIG. 4(a) which 
was digitized at a rate of 0.01 seconds and has a total length 
of 30 seconds. Therefore, the highest frequency that can be 
extracted is 25 Hz. The total cell size could be 600, but they 
have been averaged to 200 in FIG. 6(a). 

Marginal Spectrum 

The marginal spectrum offers a measure of total ampli- 
tude (or energy) contribution from each frequency value. In 
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other words, the marginal spectrum represents the cumu- 
lated amplitude over the entire data span. 

As shown in FIG. 1(a), the marginal spectrum is calcu- 
lated by the computer implemented method in step 145 after 
applying the Hilbert Transform in step 140. The marginal 5 
spectrum is the Hilbert Spectrum integrated through all time. 

In this simplification, the time coordinate is lost as in the 
Fourier spectral analysis, which leaves a summary of the 
frequency content of the event. Therefore, this presentation 
should only be used when the phenomena being analyzed is 1Q 
stationary. Formally, the marginal spectrum h(co) is defined 
as: 

r T (16) 

h(a )) = I H(o), t)d t. 

Jo 15 

Because there is no analytic expression for H(w.t), the 
integration can only be performed in a computer as a sum. 

An example of a marginal spectrum is shown in FIG. 7. 
More particularly, FIG. 7 shows the corresponding marginal 20 
spectrum of the Hilbert Spectrum given in FIG. 6(a). 

The frequency in either H(co, t) or h(co) has a totally 
different meaning from results generated by applying Fou- 
rier spectral analysis. In the Fourier representation, the 
existence of energy at a frequency, 00, means a component of 25 
a sine or a cosine wave persisted through the time span of the 
data. 

In contrast, the existence of energy at the frequency, 00 , 
means only that, in the whole time span of the data, there is 
a higher likelihood for such a wave to have appeared locally. 30 
In fact, the Hilbert Spectrum is a weighted non-normalized 
joint amplitude -frequency- time distribution. The weight 
assigned to each time -frequency cell is the local amplitude. 
Consequently, the frequency in the marginal spectrum indi- 
cates only the likelihood of an oscillation with such a 35 
frequency exists. The exact occurrence time of that oscilla- 
tion is given in the full Hilbert spectrum. 

To illustrate this difference, the corresponding Fourier 
Spectrum of the wind speed signal is also given in FIG. 7 
using a dotted line. As can be observed, there is little 40 
similarity between the Fourier spectrum and the marginal 
spectrum. While the Fourier spectrum is dominated by the 
DC term because of the non-zero mean wind speed, the 
marginal spectrum gives a nearly continuous distribution of 
energy. The Fourier spectrum is meaningless physically, 45 
because the data is not stationary. In contrast, the marginal 
spectrum provides a physically meaningful interpretation of 
the wind speed signal. 

Instantaneous Frequency 50 

There are two types of frequency modulation: the inter- 
wave and the intra-wave modulations. The first type is 
familiar because the frequency of the oscillation gradually 
changes as the waves disperse. Technically, in dispersive 
waves, the frequency is also changing within one wave, but 55 
that is generally not emphasized either for convenience, or 
for lack of a more precise frequency definition. The second 
type is less familiar, but it is also a common phenomenon: 
if the frequency changes from time to time within a wave, 
its profile can no longer be a simple sine or cosine function. 60 
Therefore, any wave profile deformation from the simple 
sinusoidal form implies intra-wave frequency modulation. 

In the past, such phenomena were treated as harmonic 
distortions. Nevertheless, such deformations should be 
viewed as intra-wave frequency modulation because the 65 
intra-wave frequency modulation is a more physically mean- 
ingful term. 
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In order to understand these frequency modulations, the 
invention applies a unique definition of instantaneous fre- 
quency. This definition stems from the EMD method and 
requires the signal to be reduced into IMF components. After 
extracting the IMF components, an instantaneous frequency 
value can be assigned to each IMF component. 
Consequently, for complicated data in which more than one 
IMF may be extracted by EMD, there can be more than one 
instantaneous frequency at a time locally. 

With the Hilbert Transform, a unique definition of instan- 
taneous frequency may be applied by the computer imple- 
mented method as illustrated by step 160. Step 160 calcu- 
lates the instantaneous frequency which is formally defined 
as follows: 

d0(t) (17) 


By calculating instantaneous frequency, step 160 of the 
invention permits the frequency value to be defined for 
every point with the value changing from point to point. 

The validity and the implications of the instantaneous 
frequency for nonlinear signals may be analyzed by exam- 
ining an idealized Stokes wave in deep water. The wave 
profile of such a Stokeian wave is modeled by 

X(f)=cos(cof+e sin c of) (18) 

Therefore, it is a intra-wave frequency modulated signal. 
Approximately, equation (18) can be shown to be: 

X(f)=(l+e/2)cos cof+e cos 2cof+ ... (19) 

The wave profile is also shown in FIG. 9(a). Because the 
intra-wave frequency can only be approximated by harmon- 
ics in Fourier analysis, we can still have the same profile, but 
not the same frequency content. The wave form shows 
sharpened crests and rounded off troughs, which make the 
crests and troughs asymmetric with respect to the mean 
surface. 

Processed with computer implemented EMD, this data 
yields only one IMF component as shown in FIG. 9(b), with 
a constant offset component (not shown). Although this 
wave has only one characteristic scale or IMF, the Wavelet 
analysis result shown in FIG. 9(c). FIG. 9(c) has many 
harmonics with two visible bands of energy corresponding 
to the highest order of approximations of the possible 
harmonics. 

In contrast, the IMF data can be processed by the inven- 
tive method to give the Hilbert Spectrum shown in FIG. 
9(d). The Hilbert Spectrum has only one frequency band 
centered around 0.03 Hz, the fundamental frequency of the 
wave train, but there is an intra-wave frequency modulation 
with a magnitude range of 0.02 to 0.04 Hz as the model truly 
represents. This intra-wave frequency modulation has been 
totally ignored in the past, for the traditional definition of 
frequency is based on the reciprocal of periodicity and 
Fourier Analysis. 

Instantaneous Energy Density 

Furthermore, the computer implemented method may also 
calculate the instantaneous energy density in step 150. The 
instantaneous energy density, IE, is defined as: 
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J' H 2 (cl>, t)du). 


( 20 ) 


Still further, this instantaneous energy density depends on 
time. Therefore, the IE can be used to check energy fluc- 
tuations. 

Stationarity 

To quantitatively measure the stationarity of a physical 10 
signal, the invention utilizes step 165 to calculate various 
stationarity measurements. Before introducing the preferred 
stationarity measurements, a brief review of conventional 
stationarity measurements is presented. 

The classic definitions of stationarity are dichotomous: a 15 
process is either stationary or nonstationary. This crude 
definition is only qualitative in nature. Such definitions are 
both overly stringent and useless at the same time: few data 
sets can satisfy the rigor of these definitions; consequently, 
no one even bothers using them to check stationarity of the 20 
signal. As a result, data as nonstationary as earthquake and 
seismological signals are routinely treated as stationary (see, 
for example, Hu, et al., 1996., Earthquake Engineering. 
Chapman & Hall, London). 

Sometimes, for some obviously nonstationary data, two 25 
less stringent definitions are invoked: piece-wise stationary 
and asymptotically stationary. These definitions are still 
dichotomous. 

To quantify the statistical processes further, an index is 
needed to give a quantitative measure of how far the process 
deviates from stationarity. A prerequisite for such a defini- 
tion is a method to present the data in the frequency-time 
space. 

With the energy-frequency-time distribution (Hilbert 35 
Spectrum) described above, stationarity of the physical 
signal may be quantitatively determined. Therefore, the 
invention introduces an index of stationarity as follows and 
calculates a Degree of Stationarity in step 165. 

The first step in defining the Degree of Stationarity, 4Q 
DS(co), is to find the mean marginal spectrum, n(oo), as 

1 ( 21 ) 

45 

Then, the Degree of Stationarity may be defined as: 


DS(a >) 




H(o), t ) 


dt. 


(22) 
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Again, the value of DS(co) can be determined by the com- 
puter. Therefore, the specialized computer 410 according to 
the invention can be treated as a stationary meter. 

For a stationary process, the Hilbert spectrum cannot be 55 
a function of time. Then, the Hilbert Spectrum will consist 
of only horizontal contour lines and DS(co) will then be 
identically zero. Only under this condition will the marginal 
spectrum be identical to the Fourier spectrum, then the 
Fourier spectrum will also make physical sense. On the other 60 
hand, if the Hilbert spectrum depends on time, the index will 
not be zero, then the Fourier spectrum will cease to make 
physical sense. 

In general, the higher the index value, the more nonsta- 
tionary is the process. The DS for the wind data is shown in 65 
FIG. 8(a). As the index shows, the data are highly nonsta- 
tionary especially for the high frequency components. 
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Eq. (22) defines the stationarity as a function of frequency. 
This is necessary because certain frequency components can 
be nonstationary while other components remain stationary. 
An example is sporadic riding wind-generated waves on 
otherwise uniform swell: the low frequency swell compo- 
nent is stationary, while the high frequency wind waves are 
intermittent, and hence nonstationary. 

The degree of stationarity can also be a function of time 
implicitly, for the definition depends on the time length of 
integration in Eq. (22). Therefore, a process can be piece- 
wise stationary. On the other hand, for a singular outburst in 
an otherwise stationary signal, the process can be regarded 
as almost stationary if a long time integral is performed, but 
nonstationary the integral only encompasses the immediate 
neighborhood of the outburst. 

Stationarity can be a complicated property of the process: 
for any signal shorter than a typical long wave period, the 
process may look transient. Yet as the signal length gets 
longer, the process can have many longer wave periods and 
becomes stationary. On the other hand, the signal can be 
locally stationary while in a long time sense nonstationary. 
An index is therefore not only useful but also necessary to 
quantify the process and give a measure of the stationarity. 

The invention also calculates a Degree of Statistic Sta- 
tionarity in step 165. The degree of stationarity defined in 
Eq. (22) can be modified slightly to include statistically 
stationary signals, for which the Degree of Statistic 
Stationarity, DSS(co,AT), is defined as 


DSS(a>, AT) = 



(23) 


where the over-line indicates averaging over a definite but 
shorter time span, AT, than the overall time duration of the 
data, T. For periodic motions, the AT can be the period. Such 
a time scale, however, is hard to define precisely for high 
dimensional, nonstationary dynamic systems. 

Even with this difficulty, the definition for DSS could be 
more useful in characterizing random variables from natural 
phenomena. Furthermore, DSS will depend on both fre- 
quency and the averaging time span. For the wind data taken 
as an example, the DSS is given in FIG. 8(a) with AT=10, 
50, 100, and 300 time steps respectively. The results show 
that while the high frequency components are nonstationary, 
they can still be statistically stationary. Two frequency bands 
at 7 and 17 Hz are highly nonstationary as the DSS averaged 
at 100 time steps shown in FIG. 8(b). These components are 
intermittent as can be seen in the IMF components and the 
marginal spectrum. A section of the original wind data is also 
plotted in FIG. 8(c) to demonstrate that there are indeed 
prominent 7 and 17 Hz time scale oscillations. 

Display of Selected Results 

The invention displays various results of the above- 
described computer implemented method in step 190. These 
displays are extremely useful in analyzing the underlying 
physics of the physical phenomenon being studied as 
described above. Furthermore, particular examples of these 
displays and the increased understanding of the underlying 
physics which these displays permit are discussed in the 
following section. 

For example, the invention generates various Hilbert 
spectra displays in the display step 190. As mentioned 
above, both color coded maps and contour maps may be 
employed to display the Hilbert spectra in display step 190. 
In addition, the color coded maps convey information to the 
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operator in a uniquely accessible way permitting a more 
thorough and deeper understanding of the physical phenom- 
enon and may be considered as necessary to analyze some 
physical phenomena. 

The displays generated by the invention in display step 
190 are central to the invention because they allow an 
operator to analyze the underlying physics of the physical 
phenomenon being studied. 

The display step 190 outputs displays to display 450. As 
mentioned above, display 450 includes devices such as a 
cathode ray tube and a flat panel display. As an alternative, 
display step 290 may generate a hard copy output by 
utilizing printer 460 or may send the generated display to 
output device 470. 

Observational Data from Laboratory and Field 
Experiments 

This section discusses a variety of geophysical energy 
signals which illustrate the usefulness of the invention. 

Earthquake Signals 

One particular application of the invention to geophysical 
signals is the analysis of earthquakes. As mentioned above, 
conventional methods cannot reveal detailed information in 
the dispersion properties, the wave form deformation, and 
the energy-frequency distribution of earthquakes because 
the signals representing the earthquake are nonlinear and 
nonstationary. 

As explained above, both nonstationarity and nonlinearity 
can induce artificial frequency smearing and reduce the true 
energy density of the earthquake signal. 

To illustrate the advantages of the computer implemented 
method disclosed herein, the following analyzes the well- 
tested data from the El Centro earthquake. 

The original El Centro earthquake signal is shown in FIG. 
10(a). These signals, when decomposed by computer imple- 
mented EMD, produce ten IMF components as shown in 
FIGS. 10(Z?)— (A:) . Furthermore, the corresponding Hilbert 
spectrum may be determined as shown in FIG. 11(a). From 
the Hilbert spectrum, one can see the diffused energy in the 
high frequency range, while persistent energy resides along 
horizontal belts below 1 Hz. 

To assure the Hilbert spectrum is a valid representation of 
the energy-frequency-time distribution, a standard Morlet 
Wavelet analysis was applied to the same data. The result is 
given in FIG. 11(6). Direct comparison between the Hilbert 
spectrum and the Wavelet analysis is not easy, for the Hilbert 
spectrum gives too many details. 

A 15xl5-point smoothing performed in step 155 results in 
the smoothed spectrum shown in FIG. 11(c). Now the 
similarity is reassuring: they show a similar pattern of 
energy concentration in the low frequency range. But the 
difference really calls for alternatives: being Fourier based, 
the Wavelet result gives much more pronounced high fre- 
quency components than does the Hilbert spectrum, a built- 
in deficiency in the Fourier analysis. 

Other than the Hilbert spectrum, the marginal spectrum 
and the corresponding Fourier spectrum are presented in 
FIG. 12(a)-(d). Although the results are presented with 
linear scales and in terms of frequency rather than period, 
the energy-frequency distribution with further smoothing is 
consistent with the previous analyses of the data shown in 
Newmark, N. M. and E. Rosenblueth, 1971, “Fundamentals 
of Earthquake Engineering ”, Prentice-Hall, Englewood 
Cliffs. N.J. The Fourier spectrum has a much wider energy 
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distribution. Consequently, the energy content for any fre- 
quency is also much lower than the corresponding values in 
the Hilbert spectrum. Furthermore, the two spectra have 
drastically different energy distribution patterns caused by 
5 the tendency of energy smearing of the Fourier spectrum. 
While both spectra give a decreasing trend in energy density 
as a function of frequency, the Fourier spectrum shows a 
conspicuous lack of energy in the low frequency range (less 
than 1 Hz), as shown in the detailed spectra covering only 
10 0 to 5 Hz frequency range. The low frequency range critical 
to high-rise buildings is severely under-represented. 

This peakiness of energy located at a very narrow fre- 
quency range can cause resonance oscillation of buildings, 
and has been observed at Mexico City during the 1985 
15 earthquake to have caused great destruction. The oscillation 
could not be duplicated by any linear model. The new 
Hilbert spectral analysis introduced here clearly identifies 
the new oscillation mode. Such oscillation modes would 
have been totally obscured by Fourier spectral analysis. 
20 From the Fourier spectra, there are no clear cut corner 
frequency or frequencies for one to construct the envelope 
spectrum (Aki and Richards, 1980, “Quantitative 
Seismology ”, W. H. Freeman, San Francisco). There are also 
no discernible stationary segments as postulated by the time 
25 history envelope method. 

Water Wave Signals 

Another application of the computer implemented method 
is water wave propagation. As mentioned above, the fre- 
quency of water waves will downshift as they propagate due 
to the interaction of weakly nonlinear wavetrains. Known 
data analysis techniques, however, has rendered proof of this 
theory nearly impossible. 

35 With the resolution power of the present invention, 
however, it is possible to show that the wave evolution 
process is not continuous and gradual, but local and discrete: 
that the well-known frequency downshift is a cumulation of 
discrete fusion of n to (n-1) waves, also known as ‘crest 
40 pairing’ (Ramamonjiarisoa and Mollo -Christensen, 1979, 
“Modulation Characteristics of Sea Surface Waves”, J. Geo- 
phys. Res., 84, 7,769-7,775. ), and Tost crest’ (Lake and 
Yuan, 1978, “A New Model for Nonlinear Gravity Waves, 
Part I., Physical Model and Experimental Evidence”, J. 
45 Fluid Mech., 88, 33-62). Some results previously reported 
by Huang, et al., 1996, “The Mechanism for Frequency 
Downshift in Nonlinear Wave Evolution”, Adv. Appl. 
Mech., 32, 59-111, illustrate this application of the Hilbert 
spectral analysis. 

50 The data were collected in the wind-wave tank at NASA 
Wallops Flight Facility (Huang and Long, 1980, “An 
Experimental Study of Surface Elevation Probability Dis- 
tribution and Statistics of Wind Generated Waves”, J. Fluid 
Mech., 101, 179-200; and Huang, et al., 1996, “The Mecha- 
55 nism for Frequency Downshift in Nonlinear Wave 
Evolution”, Adv. Appl. Mech., 32, 59-111). Surface eleva- 
tions were recorded at eight stations along the tank for a 
mechanically generated sinusoidal wave. Huang, et al., 
1996, “The Mechanism for Frequency Downshift in Non- 
60 linear Wave Evolution”, Adv. Appl. Mech., 32, 59-111 used 
the Hilbert transform and examined the change of the phase 
of the wave trains. The waves are mechanically generated by 
a wave-maker driven at 2.5 Hz with the raw data of the test 
for all eight stations given in FIG. 15(a). 

65 After performing the Hilbert transform on the data, the 
unwrapped phase functions for all eight stations are shown 
in FIG. 15(5), which indeed suggests the gradual change of 
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the phase function. The decrease in the slope of the phase 
function indicates the decrease of the mean frequency. The 
slopes of the first four stations, however, stay almost iden- 
tical to that of the first station. The first phase function that 
shows any deviation from this initially closely clustered 
group is at station 5. To this scale, the frequency change 
seems small and gradual. 

To examine the variation in more detail, take the differ- 
ence of the phase functions of all stations with respect to the 
first station and use it as an initial condition. The difference 
is shown in FIG. 15(c), from which the changes of the phase 
functions are the result of a series of steps with the sharp 
jumps confined in very short time spans. 

To quantify these phase variations, plot the phase change 
with respect to that of the first station in phase- amplitude 
diagrams as illustrated in FIGS. 15(d)-(g) for stations 2, 4, 
5 and 6, respectively. The phase- amplitude diagram for 
station 5 shows that the phase starts to jump, which reflects 
the sudden shift of the phase function shown in FIG. 15(c). 
These phase changes are very similar to the phase disloca- 
tions. As the number of jumps increases, the phase function 
seems to become a continuous and almost smooth sloped 
line. But even for those cases, the amount of phase change 
involved in each jump, as revealed by the amplitude -phase 
diagram, is still constant for all the phase jumps: 2 jt. The 2:t 
jump means a loss of exactly one wave in the process. 

Now examine the wave elevation data in detail at the jump 
points. The time series data is expanded and displayed in 
FIG. 15(A) in a two-way comparison. By superimposing the 
data from station 5 with its phase difference from station 1, 
the jump in the phase is easily associated with a single wave 
within one wave period. Then superimpose the raw elevation 
data from station 5 on those from station 1 as also shown in 
FIG. 15(A). As is apparent, all the wave peaks line up except 
at the location of the jump, where two waves become one. 
There is a loss of one wave in each phase jump event. All the 
changes are local, abrupt, and discrete. 

The Hilbert Spectrum is shown in FIG. 15(z), in which one 
can identify very local events of abrupt change of frequency 
at the 25 and 48 seconds locations. The Hilbert spectrum 
reveals the typical intra-wave frequency modulations of the 
Stokian water waves, and the more striking local shifts of 
frequency to a lower value at the locations coinciding with 
the phase jumps shown in FIG. 15(c). Another interesting 
phenomenon shows up in the marginal spectrum given in 
FIG. 1 5(f), where there is a much wider frequency distribu- 
tion than expected caused by the intra-wave frequency 
modulations. Furthermore, there is also a slight tendency of 
energy being shifted to the subharmonics frequency. 
Whether this energy flux is related to the subharmonic 
instability studied by Longuet-Higgins, 1978, Longuet- 
Higgins, M. S., 1978, “The Instabilities of Gravity Waves of 
Finite Amplitude in Deep Water, II, Subharmonics”, Proc. 
Roy. Soc. Lond., A360, 489-505. deserves further investi- 
gation. The Fourier spectrum, given in FIG. 15(A), shows 
only sidebands and harmonics without any suggestion of the 
local nature of the wave evolution. Again, the Morlet 
Wavelet analysis from the same data, shown in FIG. 15(1), 
fails to resolve any of the local changes. 

Using the computer implemented method of this 
invention, one can examine the evolution of weakly nonlin- 
ear wave trains in detail. Although the processes have 
always been assumed to be stationary globally, this inven- 
tion reveals that the processes are actually abrupt and local. 
Such variations are locally inhomogeneous; thus, they can- 
not be analyzed adequately by Fourier analysis. 
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Furthermore, the local variations are abrupt, and discrete, 
having the typical characteristics of particles. The frequency 
downshift is a cumulated effect of fusions of two waves into 
one or three waves into two. The fusion processes are 
5 irreversible. This asymmetry is due to preferred direction of 
downward energy flow. These results also raise concerns of 
the past assumptions on gradual and continuous variations of 
the wave number and frequency of a wave train. With the 
discrete jumps and the particle like nature of the individual 
10 wave, there may be a need for a new paradigm to represent 
the abrupt and discrete changes analytically. This presents a 
new challenge for nonlinear wave analysis. 

Before leaving this section, examine the data from station 
1 to illustrate the effects of the nonlinear distortion of the 
15 water wave surface. Because the emphasis is on the nonlin- 
ear distortion, the optional smoothing step 110 may be 
applied to smooth the data by a weighted running average so 
that the small irregularity will not cause mode mixing. A 
selected section of the data (5 seconds in length) and their 
20 Hilbert and Wavelet spectra are given in FIGS. 15(m) and 
15(/?) respectively. In FIG. 15(m), one can see the intra-wave 
frequency modulation. The frequency fluctuates within a 
very narrow range around 2.5 Hz, the frequency the wave- 
maker was driven. In general, the high local frequency 
25 values line up with the peaks, and the low local frequency 
values line up with the troughs as expected of the Stokes 
waves. In fact, it has been shown by Huang, et al. (1990a), 
“Wave Spectra”, The Sea, 9, 197-237; and Huang, et al. 
(1990b), “The Probability Structure of the Ocean Surface”, 
30 The Sea, 9, 335-366, that the gravity water wave can be 
approximated well by the Stokes model. On more detailed 
examination, however, the alignment is not perfect: there is 
a slight but systematic shift of the phase toward the wave 
front, an indication of the front-back asymmetry of the wave 
35 profile. 

In Wavelet result shown in FIG. 1 5(n), one cannot see any 
intra-wave frequency modulation anymore. Instead, there 
are harmonics again line up with the wave front, those 
alignments confirm the Hilbert spectrum result: the waves 
40 are front and back asymmetry. In the Wavelet spectrum, the 
energy containing frequency range is much wider than that 
is in the Hilbert spectrum, an indication of energy leakage of 
the Wavelet analysis. This comparison illustrates that the 
nonlinear wave distortion can indeed be explained by the 
45 much more physical interpretation as the intra-wave fre- 
quency modulation. 

Another application of the computer implemented method 
is to analyze ocean waves. To illustrate this application, the 
5Q following analyzes field data of ocean waves that were 
collected by the NOAA New Tidal Gauge, located at Duck, 
N.C., at the high data rate of 1 Hz. FIG. 13(a) shows the raw 
data for 60 minutes. This is typical ocean wave data from a 
field station: random and almost statistically stationary. 

55 In the past, such ocean wave data was treated with Fourier 
analysis. In fact, the studies of the wave spectra from Fourier 
analysis have been a main subject of the wave research (see, 
for example, Huang, et al., 1990a, “Wave Spectra”, The Sea, 
9, 197-237. 

60 This data when subject to the computer implemented 
EMD method yields 8 components as shown in FIGS. 
13(A)-(/), with the last component indicating the tidal varia- 
tion. The Hilbert spectrum is given in FIG. 14(a). 

For comparison, the corresponding Wavelet spectrum is 
65 shown in FIG. 14(A). While both spectra show a energy 
concentration around 0.1 Hz, the Wavelet spectrum gives a 
much more continuous distribution in time, and much wider 
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spread of energy consisting primarily of the harmonics in the 
frequency axis. 

To examine the results in detail, a short section of the data 
covering only 5 minutes is plotted in three ways. First, the 
Hilbert and the Wavelet spectra are plotted separately with 
the corresponding surface elevation of the wave data in 
FIGS. 14(c) and 14(d), respectively. In this expanded form, 
the Hilbert spectrum is very different from the Wavelet 
spectrum. Then, the two different spectra are also plotted 
together with the Hilbert spectrum in contour lines super- 
imposed on the Wavelet spectrum in FIG. 14(c). In these 
presentations, they both show the similar locations of energy 
concentration in time and frequency axes, but the Hilbert 
spectrum gives a sharper and more refined definition of the 
energy contour in FIGS. 14(c) and 14(c). Take the data and 
the spectra near 57 minute location for example: The data 
show a packet of high amplitude low frequency waves with 
frequency increasing and amplitude decreasing both before 
and after the packet. Both of these trend are vividly por- 
trayed by the Hilbert spectrum, but they are only vaguely 
suggested in the Wavelet spectrum. These variations sug- 
gested that the ocean waves are nonstationary, a conclusion 
supported by Huang, et al., 1996, “The Mechanism for 
Frequency Downshift in Nonlinear Wave Evolution”, Adv. 
Appl. Mech., 32, 59-111. 

Finally, let us also compare the marginal spectra from 
both the Hilbert and Wavelet spectra with that of Fourier 
spectrum in FIG. 14(f), in which the values of the spectra are 
staggered for clarity: the top line is the Wavelet result, the 
middle one is the Hilbert spectrum, and the bottom is the 
Fourier spectrum. The contrast is similar to the cases dis- 
cussed in the calibration and validation section: due to the 
leakage, the Wavelet spectrum is totally devoid of detail. 
While the Hilbert and direct Fourier spectra both show rich 
frequency contents, the lack of high harmonics due to either 
the nonlinear or the nonstationary effects in the Hilbert 
spectrum suggest that the Hilbert Spectrum can portray the 
energy-frequency presentation more precisely. 

Tide and Tsunami signals 

The field wave data selected here were also collected by 
a NOAA New Tidal Gauge. This gauge, located inside 
Kahului Harbor, Maui, is capable of recording data at 1 
minute steps. FIG. 16(a) shows the raw data from the tidal 
gauge for five days from Oct. 4 to Oct. 9, 1994. On Oct. 5, 
tsunami induced waves arrived at the site and created a water 
level changes of the magnitude comparable to that of the 
tidal signal. Although the tidal data is traditionally analyzed 
with Fourier expansion, the added tsunami waves are tran- 
sient. The combination, therefore, makes the whole time 
series nonstationary. Filtering cannot remove the tsunami 
signal cleanly, for the transient data and the tide will have 
many harmonic components in the same frequency range. 

The computer implemented EMD Method yields eight 
IMF components from the data as shown in FIGS. 16(Z?)-(z). 
Because the time scales of the tide and the tsunami waves 
are so different, the 8 IMF components can be easily divided 
into two groups: the high frequency signal representing the 
tsunami induced waves, and the last three low frequency 
components representing the tide. After the IMF’s were used 
to reconstitute the two separate wave motions, the raw data, 
the tidal component, and the tsunami induced waves are 
plotted together in FIG. 17. Here, the EMD serves as a filter 
to separate the tide and the transient tsunami without any 
ambiguity. 

FIG. 18(a) gives the Hilbert spectrum for both the tsunami 
and the tides from which the arrival time and the frequency 
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change of the tsunami waves are clearly shown in this 
energy-frequency-time distribution. Besides the clear dis- 
persion properties of the tsunami waves, there are two more 
interesting new observations: First, the variations of the 
5 tsunami wave frequency are phase locked with the tidal 
cycle. Second, the tsunami waves in the Harbor lasted many 
tidal cycles, with a frequency of half a cycle per hour. Such 
a period would not fit into the limited space inside the harbor 
of Kahului. The more likely explanation is for the tsunami 
10 waves being trapped in the bay, which have a much larger 
area for the wave to propagate. 

In the above presentation of FIG. 18(a), the energy in the 
tides dominates the energy in the tsunami thereby making 
the tsunami less prominent in FIG. 18(a). By using the 
15 versatile properties of the Hilbert spectrum, the invention 
can present the two phenomena separately: FIG. 1 8(b) 
shows only the tsunami components. Without the tide, the 
tsunami becomes clear. The tide components are also plotted 
over it to show the phase lock of the tsunami and the tides. 
20 FIG. 18(c) shows only the tidal components. One can see 
that in the Hilbert presentation, no harmonics are required, 
but the frequency has to be variable. 

As a demonstration of the versatility of the Hilbert 
Spectral Analysis, the tidal components can be selected to 
25 construct the Hilbert spectrum. The result is given in FIG. 
18(c), in which both the diurnal and semi-diurnal tides are 
clearly shown. 

In the marginal spectrum of the tidal components in FIG. 
30 18(d), however, the frequency for the semi-diurnal tide is 
much more variable, especially at the time when the tsunami 
just arrived, and towards the end of the measurement period. 
This variability seen in the marginal spectrum raises an 
interesting problem: whether there is interaction between the 
35 tide and the tsunami. 

Compared to the traditional Fourier spectrum, the mar- 
ginal spectrum again shows the absence of the harmonics. 
The reason seems to be clear: the tidal waves measured at a 
coastal station should be nonlinear, for the shallow water 
40 propagation is governed by a nonlinear equation. The 
absence of the harmonics raise a question concerning the 
best way to present the tidal data. Should one use harmonic 
analysis? Or should one accept the intra-wave frequency 
modulation to represent the nonlinearity? This is another 
45 interesting question for further studies. Based on the dis- 
cussion so far, it seems that the intra-wave frequency modu- 
lation representation would be more physical. 

Altimeter Data from the Equatorial Ocean 

50 Satellite altimetry has been a powerful technique for large 
scale ocean circulation studies (Huang, et al., 1978, “Ocean 
Surface Measurement Using Elevation from GEOS-3 
Altimeter”, J. Geophys. Res., 83, 4,673^1,682; Robinson, et 
al., 1983, “A Study of the Variability of Ocean Currents in 
55 the Northwestern Atlantic Using Satellite Altimetry”, J. 
Phys. Oceanogr., 13, 565-585). Because of the importance 
of the equatorial region in determining the global climate 
pattern, the altimeter data have been used extensively to 
study the dynamics of this area (Miller, et al., 1988, “GEO- 
60 SAT Altimeter Observation of Kelvin Waves and the 
1986-1987 El Nino”, Science, 239, 52-54; and Miller, et. 
al., 1990, “Large-Scale Meridional Transport in the Tropic 
Pacific Ocean During the 1986-87 El Nino From 
GEOSAT”, J. Geophys. Res. 95, 17,905-17,919; Zheng, et 
65 al., 1994, “The Effects of Shear Flow on Propagation of 
Rossby Waves in the Equatorial Oceans”, J. Phys. 
Oceanogr., 24, 1680-1686.; and Zheng, et al., 1995, “Obser- 
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vation of Equatorially Trapped Waves in the Pacific Using 
Geosat Altimeter Data”, Deep-Sea Res., (in press)). The 
accepted view of the equatorial dynamics is the propagation 
of Kelvin waves forced by variable wind stress (Byod, 1980, 
“The Nonlinear Equatorial Kelvin Waves”, J. Phys. 5 
Oceanogr., 10, 1-11; and Zheng, et al., 1995, “Observation 
of Equatorially Trapped Waves in the Pacific Using Geosat 
Altimeter Data”, Deep-Sea Res., (in press). In this model, 
the wave propagation will leave a surface elevation signa- 
ture of the order of 10 cm, which can be measured by the 10 
satellite altimeter. Such continuous data available now are 
derived from both crossover and collinear differences 
(Miller and Cheney, 1990, “Large-Scale Meridional Trans- 
port in the Tropic Pacific Ocean During the 1986-87 El Ni 
no From GEOSAT”, J. Geophys. Res. 95, 17,905-17,919), 15 
covering the period from April 1985 to September 1989. The 
final data have a spatial resolution of 8° longitude by 1° 
latitude. A typical time series on the Equator sea surface 
elevation data at 174° E. longitude is given in FIG. 19 (a). 

Limited by the data length and complicated by the 20 
dynamics, all the past investigators have encountered prob- 
lems in processing this obviously nonstationary data. 

The data has been analyzed with two methods: the Wave- 
let analysis and the Hilbert spectral analysis. Although the 
Wavelet results, given in FIGS. 19(A) shows the fluctuation 
of energy as functions of time and frequency, the energy 
distribution as function of time and frequency are too 
diffused to yield any quantitative information. 

The computer-implemented EMD method on the same 30 
sets of data yields four IMF components and a residual trend 
as shown in FIGS. 19(c)-(g). The corresponding Hilbert 
spectrum is given in FIG. 19(h). The sharpness of the 
frequency-time resolution in the Hilbert Spectrum is obvi- 
ous. The most important properties of the wave are the 35 
strong intra-wave frequency modulations. Each wave group 
will have a frequency change over 1.2 to 9.6-cycle per year 
within a period of 50 to 200 days. When we started the 
investigation of data using Hilbert spectral analysis, this was 
the first time we encountered the strong intra-wave fre- 40 
quency modulation. The curious frequency variations 
resulted in the re-examination of the classical nonlinear 
systems, which eventually clarified the role of intra-wave 
frequency modulations. (See Huang et al., 1998, The Empiri- 
cal Mode Decomposition and The Hilbert Spectrum for 45 
Nonlinear and Nonstationary Time Series Analysis, Proc. R. 
Soc. Lond. A 454, 903-995, Great Britain. 

In the present example, one can argue that, with the ocean 
depth as a small fraction of the horizontal extent, the low 
frequency geophysical wave motions are in shallow water, 50 
then the wave motions must be governed by nonlinear 
equations. Although frequency modulation phenomena are 
known to exist in Wavelet Analysis, the strength and the 
clearness of the modulations as shown here are seldom seen 
in other natural phenomena by other time-frequency distri- 55 
bution methods. 

The dynamics of the equatorial Kelvin waves has been 
discussed extensively by Zheng et al (1995). They con- 
cluded that steady wind cannot excite sustained Kelvin 
waves, but the variable wind, especially the westerly wind 60 
burst, can cause resonant interactions of the excited Kelvin 
waves with an intrinsic frequency at (1.0±0.2)xl0 -2 cpd, or 
0.22 to 0.34 cycle per year. A prevailing band of energy can 
be easily seen in FIG. 19(A). At around 400 days after the 
starting date of the data, a combination of strong low 65 
frequency (less than 0.5 cpy) and high frequency (4 to 5 
cycles per year) events occurred, which coincides with the 


onset of the 1986-1987 El Nino. Although the Hilbert 
spectrum show some energy distribution anomaly, more data 
is needed to deduce the characteristics of such special events 
as El Nino. 


Wind Signals 

Finally, let us return to wind data collected in a laboratory 
to further illustrate applications of the computer- 
implemented Sifting Process and Hilbert spectral analysis 
procedures. In the smoothed form, the Hilbert Spectrum 
looks similar to the Wavelet analysis result, yet they are still 
different. Two prominent differences stand out: the first one 
is the absence of the higher harmonics in the Hilbert 
Spectrum compared with the Wavelet result. The second one 
is the appearance of frequency modulated ridges of energy 
in the Hilbert Spectrum. Both characteristics of the Hilbert 
Spectrum are new and significant. 

In the Hilbert spectrum, we can see the increase in energy 
density with time, especially after 10 sec following the 
initiation of the data. Most of the new energy is in the 2 to 
5 Hz range with strong continuous frequency modulation 
rather than the richness of high harmonic components in the 
Wavelet analysis. As the frequency modulation is the sig- 
nature of nonlinear mechanisms, the Hilbert spectrum 
clearly shows the nonlinear properties of the turbulent flows. 

From the IMF components, we can see that c 2 through c 6 
are of the same magnitude. Collectively, they represent the 
bulk of the turbulent energy. All these components are 
highly intermittent, as indicated by the stationarity index 
study. Even after averaging, the DSS still indicates strong 
peaks at 7 and 17 Hz, which correspond to the first two peaks 
of the marginal spectrum in FIG. 7.4. These frequency 
ranges are represented by c 2 and c 6 components. In order to 
study the energy in this frequency range in detail, a special 
Hilbert spectrum for c 2 to c 6 components is constructed as 
in FIG. 20 (a). From this special Hilbert spectrum, no overall 
change is detected, but the intermittency is apparent. 
According to Liu, S. D. and S. K. Liu, 1994, “Solitary Wave 
and Turbulence ”, Shanghai Scientific and Technological 
Education Publishing House, Shanghai pg 127, nonlinearity 
and intermittency are common properties of turbulent flows. 

Next, examine the instantaneous energy density as shown 
in FIG. 20(A). As expected, the fluctuating components of 
the wind contain relatively low energy compared with the 
mean wind, or the trend, as shown by the dotted line. If the 
energy of the trend is added to the fluctuation, the overall 
energy is quite leveled, but the density is too high to be 
shown with the fluctuating components. In FIG. 20(A), the 
overall energy density has been divided by 100 in order to 
plot the two curves together for comparison. Because the 
total energy is dominated by the mean wind, it is preferably 
omitted when constructing the Hilbert spectrum. It can 
easily overwhelm all the other components; therefore, it 
should not be included unless specially justified. 

Construction of the special Hilbert spectrum with selected 
IMF components is an equivalent of band pass filtering 
through EMD, yet it is different from the Fourier band pass 
filtering, for there is no stationarity requirement as in the 
Fourier analysis. 

ALTERNATIVE EMBODIMENTS 

As described above, the invention constructs upper and 
lower envelopes 20, 30 with a cubic spline in steps 210 and 
230, respectively. This cubic spline fitting, however, has 
both overshoot and undershoot problems. These problems 
can be alleviated by using more sophisticated spline 
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methods, such as the taut spline in which the tension of the 
spline curve can be adjusted. 

Another alternative is higher-order spline fitting. 
Although such higher-order spline fitting may be more 
accurate, it will, however, introduce more inflection points 
or extrema, and consume more computing time. Therefore, 
it is not recommended as a standard operation. Only in 
special cases, it may be tried. 

As the spline fitting procedure is time consuming, more 
efficient methods can be devised by using simple mean 
values of successive extrema instead of computing the 
envelope-mean. In this way, only one spline fitting is 
required rather than two. Although this alternative is easier 
and faster to implement, the shortcomings are more severe 
amplitude averaging effects when the neighboring extrema 
are of different magnitudes. The successive-mean method 
will have a stronger forcing to reach uniform amplitudes, in 
which the true physics associated with amplitude will be 
destroyed. Therefore, the successive -me an method should 
only be applied where the amplitudes of the physical signal 
components are constants. 

Either the envelope mean or the successive mean method, 
when applied with the requirement of absolute symmetry, 
will produce the absurd result of uniform amplitude IMF’s. 
Therefore, the criteria in the Sifting Process should be 
chosen judiciously. One should avoid too stringent a crite- 
rion that we would get uniform amplitude IMF’s. On the 
other hand, one should also avoid too loose a criterion that 
would produce components deviating too much from IMF’s. 

It is well known that the most serious problem of spline 
fitting is at the ends, where cubic splines can have wide 
swings if left unattended. As an alternative, the invention 
may utilize a method of adding characteristic waves at the 
ends of the data span. This confines the large swings 
successfully. 

The method of adding characteristic waves is not con- 
ventional. In contrast, the conventional window that is often 
applied to Fourier transform data results in loss of useful 
data. To avoid this data loss and to confine swings at the ends 
of the data span, the invention extends the data beyond the 
actual data span by adding three additional characteristic 
waves at each end of the data span. 

The characteristic waves are defined by the last wave 
within the data span at the end of the data span. In other 
words, a characteristic wave is added to each end of the data 
span having an amplitude and period matching the last wave 
within the data span. This characteristic wave is a sinusoidal 
waveform that is extended three sinusoidal wave periods 
beyond the data span at each end. This process is repeated 
at the other end of the data span. In this way, spline fitting 
at the end of the data span, which can otherwise have a wide 
swing, is confined. In other words, by adding the extra 
characteristic waves at the ends beyond the data span, the 
spline curve will be tied down so that it will not have wild 
or excessive swings that would otherwise corrupt the data 
processing and analysis that utilizes these cubic splines. 

Other than the spline fitting, the Hilbert transform may 
also have end effects. Because the first and the last points of 
the data are usually of different values, the Fourier transform 
will introduce additional components to bridge over the 
difference resulting in the well-known Gibbs phenomena. To 
eliminate it in the Fourier transform, various windows have 
been adopted (see, for example, Brigham, 1974, “The fast 
Fourier Transform ”, Prentice -Hall, Englewood Cliffs, N.J.). 

Instead of a window which will eliminate some useful 
data at the end, the invention again adds two characteristic 
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waves at either end. These waves all start from zero at the 
beginning, and end at zero at the end. Thus, the annoying 
Gibbs phenomena are greatly reduced. 

Furthermore, one must be careful when analyzing weak 
5 signals embedded in a stronger signal. Both the envelope- 
mean and successive -me an depend on the existence of 
extrema. If a weak signal is imbedded in a strong signal, the 
extrema might not even be visible to the eye, but the Sifting 
Process can still pick them up. If the weak signals are phase 
10 locked with, and occur only at the maximum slope regions 
of the strong signals, then the method will have difficulties 
in picking up the extrema. In this case, the weak signals will 
appear as intra-wave frequency modulations. 

To process a weak signal embedded in a strong signal, it 
15 is necessary to separate them. This can be done by differ- 
entiating the signal once before processing. Because differ- 
entiation is a linear operator, the operation itself will not 
create nor annihilate scales. Therefore, such an operation 
could be used without adversely affecting the results. Final 
20 results then can be obtained by integration of the compo- 
nents. 

Still further, the Hilbert transform needs over-sampled 
data to define the instantaneous frequency precisely. In 
25 Fourier analysis, the Nyquist frequency is defined by two 
points per wave, but the frequency is defined for a wave 
covering the whole data span. In the invention, the instan- 
taneous frequency is defined through a differentiation 
process, and thus more data points will help defining the 
30 frequency more accurately. Based on the inventor’s 
experience, a minimum number of data points to define a 
frequency is five (or four At’s ). The lack of fine time steps 
can be alleviated through interpolating more points between 
available data. As a spline interpretation would also not 
35 create nor annihilate scales, it can also be used for the 
interpolation when the data is very jagged from under- 
sampled data. The smoothed data though have a longer 
length and are sometimes easier to process. The interpola- 
tion may give better frequency definition. 

40 Particular Limitations of The Invention 

The dependence on the existence of scale for mode 
definition has one limitation: the decomposition method 
cannot separate signals when their frequencies are too close. 
45 In this case, there would not be any characteristic scale: 
therefore, physically they are identical. This may be the most 
severe limitation of the invention, but even here, the 
invented method can still work as well as the Fourier 
Analysis. 

Particular Advantages of The Invention 

The strength of the EMD method should be reiterated. 
EMD is built on the idea of identifying the various scales in 
the data which are quantities of great physical significance. 
55 Therefore, in the Sifting Process, orthogonality is not a 
consideration, but scales are. Since orthogonal decomposi- 
tion is a characteristic for linear systems, violating this 
restriction is not a shortcoming but a breakthrough. 
Therefore, the decomposed IMF’s may or may not be 
60 orthogonal. As such, this method can be applied to nonlinear 
data. Though the IMF’s in most cases are practically 
orthogonal, it is a coincidence rather than a requirement of 
the EMD. 

Another advantage of the method is the effective use of all 
65 the data representing the physical phenomenon. In the 
Sifting Process, the longest scale is defined by the full length 
of the data. As a result, EMD can define many long period 
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oscillations. As is well known, the Hilbert transform without 
sifting tends to identify the highest frequency (Boashash, 
1992, “Estimating and Interpreting the Instantaneous Fre- 
quency of a Signal, Part I: Fundamentals”, Proc. IEEE, 80, 
520-538. ), the extraction of the long period components is 5 
indeed a new feature of the EMD. 

Finally, though the EMD method will give IMF 
components, the individual component does not guarantee 
well-defined physical meaning. This is true for all 
decompositions, especially for the methods with a priori io 
basis. In most cases, however, the IMF’s do carry physical 
significance. Great caution should be exercised in making 
such attempts. The rule for interpreting the physical signifi- 
cance of the IMF’s is that the scales should be clearly 
separated. Together with the Hilbert spectrum, the totality of 15 
the presentation should give a much more detailed repre- 
sentation of the physical processes than conventional meth- 
ods. 

The invention being thus described, it will be obvious that 
the same may be varied in many ways. Such variations are 20 
not to be regarded as a departure from the spirit and scope 
of the invention, and all such modifications as would be 
obvious to one skilled in the art are intended to be included 
within the scope of the following claims. 

I claim: 25 

1. A computer implemented method of analyzing a physi- 
cal signal representative of a physical phenomenon, com- 
prising the computer implemented steps of: 

inputting the physical signal representative of the physical 
phenomenon; 

recursively sifting the physical signal via Empirical Mode 
Decomposition to extract an intrinsic mode function 
indicative of an intrinsic oscillatory mode in the physi- 
cal phenomenon; and 35 

displaying the intrinsic mode function. 

2. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 1, further comprising the step of: 

transforming the intrinsic mode function with a Hilbert 40 
transform. 

3. The computer implemented method of analyzing a 

physical signal representative of a physical phenomenon 
according to claim 1, said recursive sifting step including the 
substeps of: 45 

identifying local maximum values in the physical signal, 
constructing an upper envelope of the physical signal 
from the identified local maximum values, 
identifying local minimum values in the physical signal, 
constructing a lower envelope of the physical signal from 
the identified local minimum values, 
determining an envelope mean from the upper and lower 
envelopes, 

generating a component signal by subtracting the enve- 55 
lope mean from the physical signal, 
treating the component signal as the physical signal, and 
recursively performing said sifting step until successive 
component signals are substantially equal. 

4. The computer implemented method of analyzing a 60 
physical signal representative of a physical phenomenon 
according to claim 3, said recursively performing substep 
including the substep of testing the component signal against 

a definition of intrinsic mode functions, 

said sifting step being recursively performed until three 65 
successive component signals satisfy the definition of 
intrinsic mode functions. 
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5. The computer implemented method of analyzing a 
physical phenomenon according to claim 3, said recursively 
performing substep including the substep of computing a 
standard deviation between successive component functions 
and comparing the standard deviation to a predetermined 
threshold value, 

said sifting step being recursively performed until the 
standard deviation falls below the predetermined 
threshold value. 

6. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 3, said recursive sifting step further 
including the substeps of: 

testing the local maximum values for an intermittency in 
the physical signal; 

said constructing an upper envelope step treating local 
maximum values failing said testing step as local 
minimum values to construct the upper envelope of the 
physical signal; 

said testing step testing the local minimum values for an 
intermittency in the physical signal; 
said constructing a lower envelope step treating local 
minimum values failing said testing step as local maxi- 
mum values to construct the lower envelope of the 
physical signal. 

7. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 1, further comprising the steps of: 

generating a residual signal by subtracting the intrinsic 
mode function from the physical signal; 
treating the residual signal as the physical signal during a 
next iteration of said recursive sifting step; 
iterating said recursive sifting step to generate a second 
intrinsic mode function indicative of a second intrinsic 
oscillatory mode in the physical phenomenon; and 
displaying the second intrinsic mode function. 

8. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 7, further comprising the steps of: 

continuing to perform said iterating step to generate an 
n-th intrinsic mode function indicative of an n-th intrin- 
sic oscillatory mode in the physical phenomenon until 
the residual signal has less than two extrema, and 
displaying the intrinsic mode functions. 

9. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 7, further comprising the step of: 

continuing to perform said iterating step to generate an 
n-th intrinsic mode function indicative of an n-th intrin- 
sic oscillatory mode in the physical phenomenon until 
the residual signal is monotonically increasing or 
decreasing, and 

displaying the intrinsic mode functions. 

10. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 8, further comprising the step of: 

applying a Hilbert transform to the intrinsic mode func- 
tion to generate a Hilbert spectrum, and 
displaying the Hilbert spectrum. 

11. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 10, further comprising the step of: 

calculating a marginal spectrum from the Hilbert 
spectrum, and 
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displaying the marginal spectrum. 

12. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 10, further comprising the step of: 

calculating an instantaneous frequency from the trans- 5 
formed intrinsic mode functions, and 
displaying the instantaneous frequency. 

13. The computer implemented method of analyzing a 

physical signal representative of a physical phenomenon 
according to claim 10, further comprising the step of: 10 

calculating an instantaneous energy density from the 
transformed intrinsic mode functions, and 
displaying the instantaneous energy density. 

14. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 15 
according to claim 10, further comprising the step of: 

calculating a measure of stationarity from the transformed 
intrinsic mode functions, and 
displaying the measure of stationarity. 

15. The computer implemented method of analyzing a 20 
physical signal representative of a physical phenomenon 
according to claim 1, further comprising the steps of: 

detecting the physical phenomenon with a sensor to 
generate an analog physical signal, and ^ 

converting the analog physical signal to a digital physical 
signal representative of the physical phenomenon. 

16. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 1, wherein the physical phenomenon is a 3Q 
geophysical phenomenon. 

17. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 16, 

wherein the geophysical phenomenon is an earthquake, an 35 
ocean wave, an ocean altitude or a wind, 
wherein the physical signal is a geophysical signal, and 
wherein said detecting step detects the geophysical 
phenomenon with a seismometer to generate a geo- 
physical signal representative of earthquake activity, a 40 
tidal gauge to generate a geophysical signal represen- 
tative of ocean wave elevation and frequency, a satellite 
altimeter to generate a geophysical signal representa- 
tive of ocean surface altitude or a pitot tube sensor to 
generate a geophysical signal representative of wind 45 
speed. 

18. The computer implemented method of analyzing a 

physical signal representative of a physical phenomenon 
according to claim 1, wherein the physical signal is nonlin- 
ear. 50 

19. The computer implemented method of analyzing a 
physical signal representative of a physical phenomenon 
according to claim 1, wherein the physical signal is nonsta- 
tionary. 

20. The computer implemented method of analyzing a 55 
physical signal representative of a physical phenomenon 
according to claim 1, wherein the physical signal is nonlin- 
ear and nonstationary. 

21. An apparatus for analyzing a physical signal repre- 
sentative of a physical phenomenon, comprising: 60 

an input device inputting the physical signal representa- 
tive of the physical phenomenon; 
a sifter recursively performing a Sifting Process on the 
physical signal using Empirical Mode Decomposition 
to extract an intrinsic mode function indicative of an 65 
intrinsic oscillatory mode in the physical phenomenon; 
and 
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a display displaying the intrinsic mode function. 

22. The apparatus for analyzing a physical signal repre- 
sentative of a physical phenomenon according to claim 21, 
further comprising: 

a Hilbert transformer transforming the intrinsic mode 
function with a Hilbert transform. 

23. The apparatus for analyzing a physical signal repre- 
sentative of a physical phenomenon according to claim 21, 
said sifter including: 

a local maximum identifier identifying local maximum 
values in the physical signal, 
an upper envelope constructor constructing an upper 
envelope of the physical signal from the identified local 
maximum values, 

a local minimum identifier identifying local minimum 
values in the physical signal, 
a lower envelope constructor constructing a lower enve- 
lope of the physical signal from the identified local 
minimum values, 

an envelope mean determiner determining an envelope 
mean from the upper and lower envelopes, 
a component signal generator generating a component 
signal by subtracting the envelope mean from the 
physical signal, 

wherein the component signal is treated as the physical 
signal during said sifter’s next recursive sifting 
Process, and 

wherein said sifter recursively performs the Sifting Pro- 
cess until successive component signals are substan- 
tially equal. 

24. The apparatus for analyzing a physical signal repre- 
sentative of a physical phenomenon according to claim 23, 
further comprising: 

a comparator comparing the component signal against a 
definition of intrinsic mode functions, 
said sifter recursively performing the Sifting Process until 
said comparator determines that three successive com- 
ponent signals satisfy the definition of intrinsic mode 
functions. 

25. The apparatus for analyzing a physical phenomenon 
according to claim 23, further comprising: 

a standard deviation calculator calculating a standard 
deviation between successive component functions, 
a comparator comparing the standard deviation to a 
predetermined threshold, 

said sifter recursively performing the Sifting Process until 
the standard deviation falls below the predetermined 
threshold value. 

26. The apparatus for analyzing a physical signal repre- 
sentative of a physical phenomenon according to claim 23, 

a tester testing the local maximum values for an intermit - 
tency in the physical signal; 

said upper envelope constructor constructing an upper 
envelope step by treating local maximum values failing 
said tester as local minimum values to construct the 
upper envelope of the physical signal; 
said tester testing the local minimum values for an inter- 
mittency in the physical signal; 
said lower envelope constructor constructing a lower 
envelope by treating local minimum values failing said 
tester as local maximum values to construct the lower 
envelope of the physical signal. 

27. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 21, further com- 
prising: 
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a residual signal generator generating a residual signal by 
subtracting the intrinsic mode function from the physi- 
cal signal; 

wherein the residual signal is treated as the physical signal 
during a next iteration of the Sifting Process performed 5 
by said sifter; 

an iterator iterating the Sifting Process performed by said 
sifter to generate a second intrinsic mode function 
indicative of a second intrinsic oscillatory mode in the 
physical phenomenon; and 10 

a display displaying the second intrinsic mode function. 

28. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 27, wherein said 
sifter continues to perform said Sifting Process to generate 

an n-th intrinsic mode function indicative of an n-th intrinsic ^ 
oscillatory mode in the physical phenomenon until the 
residual signal has less than two local extrema, and 

wherein said display displays the intrinsic mode func- 
tions. 

29. The apparatus for analyzing physical data from a 20 
physical phenomenon according to claim 27, wherein said 
sifter continues to perform said Sifting Process to generate 

an n-th intrinsic mode function indicative of an n-th intrinsic 
oscillatory mode in the physical phenomenon until the 
residual signal is monotonically increasing or decreasing, 25 
and 

wherein said display displays the intrinsic mode func- 
tions. 

30. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 28, further com- 30 
prising: 

a Hilbert Spectrum generator applying a Hilbert transform 
to the intrinsic mode function to generate a Hilbert 
spectrum, 

wherein said display displays the Hilbert spectrum. 35 

31. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 30, further com- 
prising: 

a marginal spectrum calculator calculating a marginal 
spectrum from the Hilbert spectrum, 40 

wherein said display displays the marginal spectrum. 

32. The apparatus for analyzing physical data from a 

physical phenomenon according to claim 30, further com- 
prising: 45 

an instantaneous frequency calculator calculating an 
instantaneous frequency from transformed intrinsic 
mode functions, 

wherein said display displays the instantaneous fre- 
quency. 50 

33. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 30, further com- 
prising: 

an instantaneous energy density calculator calculating an 
instantaneous energy density from the transformed 55 
intrinsic mode functions, 

wherein said display displays the instantaneous energy 
density. 

34. The apparatus for analyzing physical data from a 

physical phenomenon according to claim 30, further com- 60 

prising: 

a stationarity calculator calculating a measure of station- 
arity from the transformed intrinsic mode functions, 
wherein said display displays the measure of stationarity. 

35. The apparatus for analyzing physical data from a 65 
physical phenomenon according to claim 21, further com- 
prising: 
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a sensor detecting the physical phenomenon to generate 
an analog physical signal, and 
an analog to digital convertor converting the analog 
physical signal to a digital physical signal representa- 
tive of the physical phenomenon. 

36. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 21, wherein the 
physical phenomenon is a geophysical phenomenon. 

37. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 36, wherein the 
geophysical phenomenon is an earthquake, an ocean wave, 
an ocean altitude or a wind and wherein the physical signal 
is a geophysical signal reproduction of the geophysical 
phenomenon, the apparatus further comprising: 

a seismometer generating a geophysical signal represen- 
tative of earthquake activity, a tidal gauge generating a 
geophysical signal representative of ocean wave eleva- 
tion and frequency, a satellite altimeter generating a 
geophysical signal representative of ocean surface 
altitude, or a pitot tube sensor generating a geophysical 
signal representative of wind speed. 

38. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 21, wherein the 
physical signal is nonlinear. 

39. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 21, wherein the 
physical signal is nonstationary. 

40. The apparatus for analyzing physical data from a 
physical phenomenon according to claim 21, wherein the 
physical signal is nonlinear and nonstationary. 

41. An article of manufacture, comprising: 

a computer-usable medium including computer-readable 
program code means, embodied therein, for causing a 
computer to analyze a physical signal representative of 
a physical phenomenon, the computer-readable pro- 
gram code means comprising: 
computer-readable program code means for inputting the 
physical signal representative of the physical phenom- 
enon; 

computer-readable program code means for recursively 
sifting the physical signal via Empirical Mode Decom- 
position to extract an intrinsic mode function indicative 
of an intrinsic oscillatory mode in the physical phe- 
nomenon; and 

computer-readable program code means for displaying 
the intrinsic mode function on a display. 

42. The article of manufacture according to claim 41, the 
computer-readable program code means further comprising: 

computer-readable program code means for transforming 
the intrinsic mode function with a Hilbert transform. 

43. The article of manufacture according to claim 41, said 
recursive sifting means including: 

computer-readable program code means for identifying 
local maximum values in the physical signal, 
computer-readable program code means for constructing 
an upper envelope of the physical signal from the 
identified local maximum values, 
computer-readable program code means for identifying 
local minimum values in the physical signal, 
computer-readable program code means for constructing 
a lower envelope of the physical signal from the 
identified local minimum values, 
computer-readable program code means for determining 
an envelope mean from the upper and lower envelopes, 
computer-readable program code means for generating a 
component signal by subtracting the envelope mean 
from the physical signal, 
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computer-readable program code means for treating the 
component signal as the physical signal, and 

computer-readable program code means for recursively 
operating said sifting means until successive compo- 
nent signals are substantially equal. 

44. The article of manufacture according to claim 43, said 
recursively operating means including: 

computer-readable program code means for comparing 
the component signal to determine whether the com- 
ponent signal satisfies a definition of Intrinsic Mode 
Functions, 

wherein the operation performed by said sifting means are 
recursively performed until three successive compo- 
nent signals satisfy the definition of Intrinsic Mode 
Functions. 

45. The article of manufacture according to claim 43, said 
recursively operating means including: 

computer-readable program code means for computing a 
standard deviation between successive component 
functions and comparing the standard deviation to a 
predetermined threshold value, 

wherein the operations performed by said sifting means 
are recursively performed until the standard deviation 
falls below the predetermined threshold value. 

46. The article of manufacture according to claim 43, said 
recursive sifting means further including: 

computer-readable program code means for testing the 
local maximum values for an intermittency in the 
physical signal; 

said constructing an upper envelope means treating local 
maximum values failing said testing step as local 
minimum values to construct the upper envelope of the 
physical signal; 

said testing means testing the local minimum values for 
an intermittency in the physical signal; 

said constructing a lower envelope means treating local 
minimum values failing said testing step as local maxi- 
mum values to construct the lower envelope of the 
physical signal. 

47. The article of manufacture according to claim 41, the 
computer-readable program code means further comprising: 

computer-readable program code means for generating a 
residual signal by subtracting the intrinsic mode func- 
tion from the physical signal; 

computer-readable program code means for treating the 
residual signal as the physical signal during a next 
iteration of the operations performed by said recursive 
sifting means; 

computer-readable program code means for iterating the 
operations performed by said recursive sifting means to 
generate a second intrinsic mode function indicative of 
a second intrinsic oscillatory mode in the physical 
phenomenon; and 

computer-readable program code means for displaying 
the second intrinsic mode function. 

48. The article of manufacture according to claim 47, the 
computer-readable program code means further comprising: 

computer-readable program code means for continuing 
the operations performed by said iterating means to 
generate an n-th intrinsic mode function indicative of 
an n-th intrinsic oscillatory mode in the physical phe- 
nomenon until the residual signal contains less than two 
local extrema, and 
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computer-readable program code means for displaying 
the intrinsic mode functions. 

49. The article of manufacture according to claim 47, the 
computer-readable program code means further comprising: 

5 computer-readable program code means for continuing 
the operations performed by said iterating means to 
generate an n-th intrinsic mode function indicative of 
an n-th intrinsic oscillatory mode in the physical phe- 
nomenon until the residual signal is monotonically 
increasing or decreasing, and 

computer-readable program code means for displaying 
the intrinsic mode functions. 

50. The article of manufacture according to claim 48, the 
15 computer-readable program code means further comprising: 

computer-readable program code means for applying a 
Hilbert transform to the intrinsic mode function to 
generate a Hilbert spectrum, and 
20 computer-readable program code means for displaying 
the Hilbert spectrum. 

51. The article of manufacture according to claim 50, the 
computer-readable program code means further comprising: 

computer-readable program code means for calculating a 
25 marginal spectrum from the Hilbert spectrum, and 

computer-readable program code means for displaying 
the marginal spectrum. 

52. The article of manufacture according to claim 50, the 
computer-readable program code means further comprising: 

computer-readable program code means for calculating an 
instantaneous frequency from the transformed intrinsic 
mode functions, and 

computer-readable program code means for displaying 
35 the instantaneous frequency. 

53. The article of manufacture according to claim 50, the 
computer-readable program code means further comprising: 

computer-readable program code means for calculating an 
4Q instantaneous energy density from the transformed 

intrinsic mode functions, and 

computer-readable program code means for displaying 
the instantaneous energy density. 

54. The article of manufacture according to claim 50, the 
45 computer-readable program code means further comprising: 

computer-readable program code means for calculating a 
measure of stationarity from the transformed intrinsic 
mode functions, and 

50 computer-readable program code means for displaying 
the measure of stationarity. 

55. The article of manufacture according to claim 41, 
wherein the physical signal is a geophysical signal repre- 
sentative of a geophysical phenomenon. 

55 56. The article of manufacture according to claim 55 

wherein the geophysical phenomenon is an earthquake, an 
ocean wave, an ocean altitude, or a wind. 

57. The article of manufacture according to claim 41, 
wherein the physical data from physical signal is nonliner. 
60 58. The article of manufacture according to claim 41, 

wherein the physical data from physical signal is nonsta- 
tionary. 

59. The article of manufacture according to claim 41, 
wherein the physical signal is nonlinear and nonstationary. 





